CN111507415B - 一种基于分布密度的多源大气数据聚类方法 - Google Patents

一种基于分布密度的多源大气数据聚类方法 Download PDF

Info

Publication number
CN111507415B
CN111507415B CN202010314605.3A CN202010314605A CN111507415B CN 111507415 B CN111507415 B CN 111507415B CN 202010314605 A CN202010314605 A CN 202010314605A CN 111507415 B CN111507415 B CN 111507415B
Authority
CN
China
Prior art keywords
matrix
cid
row
vertex
edge
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010314605.3A
Other languages
English (en)
Other versions
CN111507415A (zh
Inventor
樊仲欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Information Science and Technology
Original Assignee
Nanjing University of Information Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Information Science and Technology filed Critical Nanjing University of Information Science and Technology
Priority to CN202010314605.3A priority Critical patent/CN111507415B/zh
Publication of CN111507415A publication Critical patent/CN111507415A/zh
Application granted granted Critical
Publication of CN111507415B publication Critical patent/CN111507415B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/23Clustering techniques
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/90Details of database functions independent of the retrieved data types
    • G06F16/901Indexing; Data structures therefor; Storage structures
    • G06F16/9024Graphs; Linked lists
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/90Details of database functions independent of the retrieved data types
    • G06F16/906Clustering; Classification
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Software Systems (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computing Systems (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Algebra (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Complex Calculations (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

本发明公开了一种基于分布密度的多源大气数据聚类方法,首先,构建一个M维数据组成的数据量为N的数据集DS,并判断数据集DS的聚类趋势;其次,生成数据集的距离矩阵DM的全邻域分布密度矩阵DDM;然后,以分布密度阈值ddth为参数,划分出全邻域分布密度矩阵DDM的密度峰值与离散点;最后,截取全体数据的边矩阵E将部分离散点合并入密度峰值,得到聚类结果。本发明仅使用分布密度阈值这一单一参数实现控制聚类结果,并确保可聚类任意分布形状、分布均匀性的数据;且能自动分离出噪点。

Description

一种基于分布密度的多源大气数据聚类方法
技术领域
本发明属于数据挖掘领域,具体涉及一种基于分布密度的多源大气数据聚类方法。
背景技术
在大数据挖掘和分析的实际应用中,数据是从不同领域的不同来源收集或从不同的特征采集器获取的。例如,网站上共享的某张图片往往有不同来源的文本标记和描述;特定的新闻由多个新闻机构报道;相同的语义(如hello)用多种语言形式表示;图像由不同类型的特征描述。所有这些都称为多源数据(或多视图数据)。这些数据表现出异质性,但又具有潜在的联系。换言之,在这些数据中每个单独的源(或视图)对于知识发现任务都有其特定的属性,而不同的源通常包含应加以利用的互补信息。因此,如何利用这些信息挖掘多源数据的潜在价值,在大数据研究中具有十分重要的意义。在现实数据分析中的应用也需要能够处理多源数据对象的先进技术,以将数据挖掘和知识发现推向新的高度。
大气数据是具有多源性的典型范例,这是因为近现代科学技术的迅猛发展导致了对于大气的观测和监测手段不断趋于多样化,出现了地基(自动站、雷达)、海基(浮标)、空基(气球、飞机)、天基(卫星)这样一套全方位的立体化观测监测体系,而不同的观测监测手段所产生的多源数据是既各有特点又相互关联的,所以对多源大气数据进行聚类分析,根据数据的属性性质进行归类,这对于类簇的发掘(常见天气现象统计)、离群点或噪点检测(灾害性天气识别)以及更进一步的天气预测预警等都是必不可少的步骤。
聚类算法主要有划分式、密度式,网格式,层次式,其中划分式聚类以K-means算法为代表,通过指定簇数以及寻找簇中心点进行聚类,密度式聚类以DBSCAN算法为代表,通过查找密度相连对象的最大集合进行聚类,网格式聚类以CLIQUE算法为代表,它把每个数据维划分成不重叠的网格,从而把数据对象整个嵌入网格空间划分成单元再进行基于密度的聚类,层次式聚类以BIRCH算法为代表,它建立一个聚类特征树,然后对树的叶子结点按簇直径进行分裂从而得到聚类结果。上述各类聚类算法经过不断的发展和相互借鉴融合后出现了一种经典的密度峰值聚类算法DPC,该算法刻画聚类中心为本身密度大并且与其它密度更大的数据点间的距离相对更大,进而实现了用单一的截断距离参数控制任意分布形状的数据高效聚类,但是并不是所有的数据集都能通过决策图准确地找到簇中心而且该算法并不能判别噪点,因此有研究人员是对DPC算法进行改进,力求解决该算法截断距离的确定和簇中心的选择这两大问题,而这虽然取得了一定的成效,但簇中心选取及未分配点分配准确率却仍然不高,一方面原因是簇中心难以准确判明,另一方面是由于其所采取的未分配点最近邻分配策略不完善,因此,本发明提出一种全新的基于数据全邻域分布密度的多源大气数据聚类分析算法,该算法可以解决现有技术中存在的增加聚类参数和无法自动化的判别噪点的问题。
发明内容
发明目的:本发明提供一种基于分布密度的多源大气数据聚类方法,可通过单一参数实现控制聚类结果,且能自动找出簇边界点并分离出噪点。
技术方案:本发明所述的一种基于分布密度的多源大气数据聚类方法,包括以下步骤:
(1)构建数据集DS,并判断数据集DS的聚类趋势;所述数据集DS为一个M维数据组成的数据量为N的数据集;
(2)生成数据集的距离矩阵DM的全邻域分布密度矩阵DDM;
(3)以分布密度阈值ddth为参数,划分出全邻域分布密度矩阵DDM的密度峰值与离散点;
(4)截取全体数据的边矩阵E将簇边界点合并入簇核,得到聚类结果。
进一步地,所述步骤(1)包括以下步骤:
(11)数据集xi=[xi1...xiM]i∈{1,...N}表示数据集的第i个顶点,idi=i表示第i个顶点的唯一标识;计算距离矩阵/>其中D(xi,xj)i,j∈{1,...N}表示两个顶点间的距离,以及距离矩阵的另一种表示形式其中i,j∈{1,...N}且i≠j表示第i个和第j个顶点的唯一标识,D(xi,xj)i,j∈{1,...N}且i≠j表示两个顶点间的距离;
(12)从DE中查找第三列的最小值,将其所在行,即边,从DE删除后加入空的边矩阵E;
(13)不断从DE中查找一个顶点唯一标识属于E同时另一顶点唯一标识属于DE顶点唯一标识集合减去E顶点唯一标识集合的所有边,再将其中第三列最小值的边从DE删除后按行拼接加入E,直到DE的顶点唯一标识集合包含于E的顶点唯一标识集合为止,得到最小距离连通的边矩阵E;
(14)生成最小距离连通图MDG[ID,E],其中ID={id1,id2,…idN},然后取边矩阵E的第三列为DECNID
(15)计算DECNID中除第一行外每一行和上一行的差值,取其中的最大值Dmax在DECNID中对应行的值和对应行下一行的值计算平均值,得到均值deth,deth即为肘阈值;
(16)用DECNID中小于deth的值建立距离集合Dsmall,用DECNID中大于deth的值建立距离集合Dbig;计算出距离集合Dsmall、Dbig中元素的个数numsmall、numbig;均值meansmall、meanbig,以及Dsmall、DECNID的变异系数cvsmall、cvECNID
(17)根据以下公式得到聚类趋势指数CTI:
CTI∈[-1,1)
式中,
(18)如CTI>0.5则进入步骤(2),否则提示“无聚类趋势”信息并结束本次聚类流程。
进一步地,所述步骤(2)包括以下步骤:
(21)以距离矩阵DM的各元素为半径r,通过超球体公式和维度M得到超球体体积矩阵VM;
(22)对体积矩阵VM的第1行VMr1按由小到大顺序排序得到VMsr1以及VMr1元素顺序在排序后的重新排列Ir1,然后按照i={1,2…N}的顺序遍历VMsr1的体积值VMsr1(i),计算累计量同时生成VMsr1(i)对应的顶点数量PNsr1(i)=i,但如果遇到VMsr1(i)=VMsr1(i-1)=…=VMsr1(i-j)的情况则PNsr1(i)=PNsr1(i-1)=…=PNsr1(i-j)=i;
(23)将向量AVMsr1,PNsr1按照VMr1的元素顺序Ir1重新排回原来顺序得到AVMr1,PNr1
(24)对体积矩阵VM的第2至N行VMr2,VMr3…VMrN执行步骤(22)至步骤(23),得到Ir2,…IrN和AVMr2,…AVMrN和PNr2,…PNrN
(25)将向量集合{Ir1,Ir2,…IrN}依序逐行组合成行元素顺序排序矩阵I,I为N×N非对称阵,将向量集合{AVMr1,AVMr2,…AVMrN}依序逐行组合成超球体累计体积矩阵AVM,AVM为N×N非对称阵,将向量集合{PNr1,PNr2,…PNrN}依序逐行组合成顶点数量矩阵PN,PN为N×N非对称阵,计算得到DDM,DDM为N×N非对称矩阵如下:
进一步地,所述步骤(3)包括以下步骤:
(31)设待标记顶点唯一标识集合AID=ID,密度峰值的顶点唯一标识集合离散点的顶点唯一标识集合/>并创建空队列QU;
(32)从AID中选择第1个顶点唯一标识AID(1)并查找其在DDM中的对应行DDMrAID(1),将DDMrAID(1)按IrAID(1)重新排序得到顺序排列的DDMsrAID(1),然后从第2个元素开始依序遍历DDMsrAID(1)找到前ci个大于阈值ddth的元素,如果ci=0则将AID(1)从AID删除再加入到NID集合中并跳转至步骤(35),否则将ci个元素在DDM中对应列的唯一标识ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))和AID取交集后入队列QU,并将AID(1),ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))从AID删除再加入到新建集合CIDnew={AID(1),ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))};
(33)如QU为空则进入步骤(34),否则顶点唯一标识ID(i)出队列QU,查找ID(i)在DDM中的对应行DDMrID(i),将DDMrID(i)按IrID(i)重新排序得到DDMsrID(i),然后从第2个元素开始依序遍历DDMsrID(i)找到前ci个大于阈值ddth的元素,如果ci=0跳转至本步骤开头,否则将ci个元素在DDM中对应列的唯一标识ID(IrID(i)(2))…ID(IrID(i)(ci+1))和AID取交集后入队列QU,并将ID(i),ID(IrID(i)(2))…ID(IrID(i)(ci+1))从AID删除再更新集合CIDnew=CIDnew∪{ID(i),ID(IrID(i)(2))…ID(IrID(i)(ci+1))},跳转至本步骤开头;
(34)检查最新生成并更新完成的集合CIDnew,如果CIDnew和CID集合中的某些子集合CIDj,…CIDk的交集不为空,则CIDnew=CIDnew∪CIDj∪…CIDk,并将CIDj,…CIDk从CID删除后将CIDnew加入到CID成为新增的子集合;
(35)如AID不为空则转到步骤(32)继续执行,否则得到密度峰值顶点集合CID和离散顶点集合NID并进入步骤(4)。
进一步地,所述步骤(4)包括以下步骤:
(41)创建N×N矩阵EM赋值-∞,然后遍历N个顶点的MDG[ID,E]中边矩阵E的各行,即各条边eqp=[idq,idp,dqp]令得到E的另一种表示形式EM;
(42)从距离矩阵DM中截取离散顶点集合NID为行序号,密度峰值顶点集合CID子集合的并集为列序号的部分,求其最小值,得到其中NID(row)表示最小值对应的离散顶点唯一标识,CID(col)表示最小值对应的密度峰值顶点唯一标识,而其所属的CID子集合表示为CIDcol
(43)从DM、EM中截取CNID=CIDcol∪NID(row)部分,得到CDM=DM(CNID,CNID)、CEM=EM(CNID,CNID);取CEM中非-∞元素组成向量DECNID,如果DECNID的元素数量小于CNID的元素数量减一,则进入步骤(44),否则进入步骤(45);
(44)将CEM因截取导致的断开部分补全,重建DECNID
(45)执行步骤(15)至(17);
(46)如CTI≤0.5则将NID(row)从NID删除再加入到CIDcol并转到步骤(42)继续执行,否则得到最终聚类结果簇集合CID包括密度峰值与簇边界点,以及噪点集合NID,结束所有步骤。
进一步地,所述步骤(44)包括以下步骤:
(441)从CEM中取任一非-∞元素d1i1j生成一条边CE1=[id1i,id1j,d1i1j],其中id1i,id1j表示d1i1j对应CEM的行序号和列序号,令EXI={id1i,id1j},创建空队列QU2并将id1i,id1j入队列;
(442)如QU2为空则进入步骤(443),否则顶点唯一标识idi出队列QU2,在CEM中查找与idi在同行和同列的排除EXI元素的列序号和行序号idj,…idk及其非-∞元素dj,…dk,组成新建矩阵然后令/>EXI=EXI∪{idj,…idk},将idj,…idk入队列QU2,将CE1i按行拼接加入CE1生成边矩阵CE1,跳转至本步骤开头;
(443)对CEM执行步骤(441)至步骤(442),每次执行生成新的边矩阵CE2,…CEn,直到CEM中所有元素均为-∞为止,得到边矩阵集合CE={CE1,CE2,…CEn};
(444)从距离矩阵CDM中查找出CE里面各断开部分之间的最小距离,得到边矩阵其中ICE的第三列表示CE里面每两个边矩阵的顶点唯一标识集合(边矩阵第一二列的并集)之间的最小距离,而第一二列分别为其对应的CE中各边矩阵的序号,即边矩阵唯一标识;
(445)从ICE中查找第三列的最小值,将其所在边从ICE删除后加入空的边矩阵CEn+1
(446)不断从ICE中查找一个边矩阵唯一标识属于CEn+1同时另一边矩阵唯一标识属于ICE边矩阵唯一标识集合减去CEn+1边矩阵唯一标识集合的所有边,再将其中第三列最小值的边从ICE删除后按行拼接加入CEn+1,直到ICE的边矩阵唯一标识集合包含于CEn+1的边矩阵唯一标识集合为止,得到CE中各边矩阵间断开部分的补全边矩阵CEn+1
(447)按行拼接边矩阵CE1,CE2,…CEn,CEn+1得到CE2=[CE1;CE2;…CEn;CEn+1],取CE2的第三列为向量DECNID
有益效果:与现有技术相比,本发明的有益效果:1、本发明仅使用分布密度阈值这一单一参数实现控制聚类结果,并确保可聚类任意分布形状、分布均匀性的数据;2、本发明能自动找出簇边界点并分离出噪点:本发明提出的ANDD算法先用簇分布密度判定的方式划分出密度峰值,即簇核;与离散点,即簇边界点和噪点,完成初步聚类以后,再用聚类趋势分析算法将簇边界点逐步合并入簇核,从而分离出噪点并确定最终的聚类结果。
附图说明
图1为基于分布密度的多源大气数据聚类模型结构示意图;
图2为ANDD(ddth=4.2)算法在eye数据集上的聚类效果图;
图3为ANDD(ddth=0.8)算法在flame数据集上的聚类效果图;
图4为KNN-DPC(K=5)算法在flame数据集上的聚类效果图;
图5为DPC(dc=5)算法在flame数据集上的聚类效果图;
图6为AP(p=35)算法在flame数据集上的聚类效果图;
图7为DBSCAN(Eps=0.09/Minpts=8)算法在flame数据集上的聚类效果图;
图8为ANDD(ddth=8×10-9)算法在s2数据集上的聚类效果图;
图9为KNN-DPC(K=12)算法在s2数据集上的聚类效果图;
图10为DPC(dc=2)算法在s2数据集上的聚类效果图;
图11为AP(p=20)算法在s2数据集上的聚类效果图;
图12为DBSCAN(Eps=0.23/Minpts=7)算法在s2数据集上的聚类效果图。
具体实施方式
下面结合附图对本发明作进一步详细说明。
分布密度为在一个M维数据组成的数据量为N的数据集DS中,任意一项数据(顶点)i到另一项数据j的分布密度dd(i,j)为:
其中V(i,k)表示以顶点i为球心,i到k的距离为半径的超球体体积,DS(i,j)表示V(i,j)范围内的顶点,PN(i,j)表示V(i,j)范围内的顶点数量,超球体公式如下:
其中r为半径,M为数据维度,Γ为伽玛函数。
上述定义分布密度的公式(1)表示的是以数据集中任意两个顶点距离为半径所形成的超球体内,顶点数量与球心到所有顶点的超球体体积总和之比,因此在该定义中顶点的体积总和越小(也就是顶点的分布离球心越接近越集中)则分布密度就越大。
本发明提出一种全新的基于数据全邻域分布密度的多源大气数据聚类分析算法(All Neighborhood Distribution Density clustering algorithm,ANDD),以分布密度dd(i,j)的阈值ddth为判别标准,对数据集DS和距离矩阵DM(以距离值作为元素的N×N对称二维数组)进行聚类的流程步骤如下:
步骤1:构建数据集DS,并判断数据集DS的聚类趋势。
(1.1)数据集其中M为数据维度,N为数据量,xi=[xi1...xiM]i∈{1,...N}表示数据集的第i个顶点,idi=i表示第i个顶点的唯一标识;计算距离矩阵/>其中D(xi,xj)i,j∈{1,...N}表示两个顶点间的距离,以及距离矩阵的另一种表示形式/>其中i,j∈{1,...N}且i≠j表示第i个和第j个顶点的唯一标识,D(xi,xj)i,j∈{1,...N}且i≠j表示两个顶点间的距离。
(1.2)从DE中查找第三列的最小值,将其所在行,即边,从DE删除后加入空的边矩阵E。
(1.3)不断从DE中查找一个顶点唯一标识属于E同时另一顶点唯一标识属于DE顶点唯一标识集合减去E顶点唯一标识集合的所有边,再将其中第三列最小值的边从DE删除后按行拼接加入E,直到DE的顶点唯一标识集合包含于E的顶点唯一标识集合为止,得到最小距离连通的边矩阵E。
(1.4)生成最小距离连通图MDG[ID,E],其中ID={id1,id2,…idN},然后取边矩阵E的第三列为DECNID
(1.5)计算DECNID中除第一行外每一行和上一行的差值,取其中的最大值Dmax在DECNID中对应行的值和对应行下一行的值计算平均值,得到均值deth,deth即为肘阈值。
(1.6)用DECNID中小于deth的值建立距离集合Dsmall,用DECNID中大于deth的值建立距离集合Dbig;计算出距离集合Dsmall、Dbig中元素的个数numsmall、numbig;均值meansmall、meanbig,以及Dsmall、DECNID的变异系数cvsmall、cvECNID
(1.7)根据以下公式得到聚类趋势指数CTI:
CTI∈[-1,1)
式中,
(1.8)如CTI>0.5则进入步骤2,否则提示“无聚类趋势”信息并结束本次聚类流程。
步骤2:生成数据集的距离矩阵DM的全邻域分布密度矩阵DDM。
(2.1)以距离矩阵DM的各元素为半径r,通过公式(2)和维度M得到超球体体积矩阵VM。
(2.2)对体积矩阵VM的第1行VMr1按由小到大顺序排序得到VMsr1以及VMr1元素顺序在排序后的重新排列Ir1,然后按照i={1,2…N}的顺序遍历VMsr1的体积值VMsr1(i),计算累计量同时生成VMsr1(i)对应的顶点数量PNsr1(i)=i,但如果遇到VMsr1(i)=VMsr1(i-1)=…=VMsr1(i-j)的情况则PNsr1(i)=PNsr1(i-1)=…=PNsr1(i-j)=i。
(2.3)将向量AVMsr1,PNsr1按照VMr1的元素顺序Ir1重新排回原来顺序得到AVMr1,PNr1
(2.4)对体积矩阵VM的第2至N行VMr2,VMr3…VMrN执行步骤(2.2)至步骤(2.3),得到Ir2,…IrN和AVMr2,…AVMrN和PNr2,…PNrN
(2.5)将向量集合{Ir1,Ir2,…IrN}依序逐行组合成行元素顺序排序矩阵I,I为N×N非对称阵,将向量集合{AVMr1,AVMr2,…AVMrN}依序逐行组合成超球体累计体积矩阵AVM,AVM为N×N非对称阵,将向量集合{PNr1,PNr2,…PNrN}依序逐行组合成顶点数量矩阵PN,PN为N×N非对称阵,计算得到DDM,DDM为N×N非对称矩阵如式(3):
步骤3:以分布密度阈值ddth为参数,划分出全邻域分布密度矩阵DDM的密度峰值,即簇核;离散点,即簇边界点和噪点。
(3.1)设待标记顶点唯一标识集合AID=ID,密度峰值的顶点唯一标识集合离散点的顶点唯一标识集合/>并创建空队列QU。
(3.2)从AID中选择第1个顶点唯一标识AID(1)并查找其在DDM中的对应行DDMrAID(1),将DDMrAID(1)按IrAID(1)重新排序得到顺序排列的DDMsrAID(1),然后从第2个元素开始依序遍历DDMsrAID(1)找到前ci个大于阈值ddth的元素,如果ci=0则将AID(1)从AID删除再加入到NID集合中并跳转至步骤(3.5),否则将ci个元素在DDM中对应列的唯一标识ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))和AID取交集后入队列QU,并将AID(1),ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))从AID删除再加入到新建集合CIDnew={AID(1),ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))}。
(3.3)如QU为空则进入步骤(3.4),否则顶点唯一标识ID(i)出队列QU,查找ID(i)在DDM中的对应行DDMrID(i),将DDMrID(i)按IrID(i)重新排序得到DDMsrID(i),然后从第2个元素开始依序遍历DDMsrID(i)找到前ci个大于阈值ddth的元素,如果ci=0跳转至本步骤开头,否则将ci个元素在DDM中对应列的唯一标识ID(IrID(i)(2))…ID(IrID(i)(ci+1))和AID取交集后入队列QU,并将ID(i),ID(IrID(i)(2))…ID(IrID(i)(ci+1))从AID删除再更新集合CIDnew=CIDnew∪{ID(i),ID(IrID(i)(2))…ID(IrID(i)(ci+1))},跳转至本步骤开头;
(3.4)检查最新生成并更新完成的集合CIDnew,如果CIDnew和CID集合中的某些子集合CIDj,…CIDk的交集不为空,则CIDnew=CIDnew∪CIDj∪…CIDk,并将CIDj,…CIDk从CID删除后将CIDnew加入到CID成为新增的子集合。
(3.5)如AID不为空则转到步骤(3.2)继续执行,否则得到密度峰值顶点集合CID和离散顶点集合NID并进入步骤4。
步骤4:截取全体数据的边矩阵E将簇边界点合并入簇核,得到聚类结果。
(4.1)创建N×N矩阵EM赋值-∞,然后遍历N个顶点的MDG[ID,E]中边矩阵E的各行,即各条边eqp=[idq,idp,dqp]令得到E的另一种表示形式EM。
(4.2)从距离矩阵DM中截取离散顶点集合NID为行序号,密度峰值顶点集合CID子集合的并集为列序号的部分,求其最小值,得到其中NID(row)表示最小值对应的离散顶点唯一标识,CID(col)表示最小值对应的密度峰值顶点唯一标识而其所属的CID子集合表示为CIDcol
(4.3)从DM、EM中截取CNID=CIDcol∪NID(row)部分,得到CDM=DM(CNID,CNID)、CEM=EM(CNID,CNID)。取CEM中非-∞元素组成向量DECNID,如果DECNID的元素数量小于CNID的元素数量减一,则进入步骤(4.4),否则进入步骤(4.5)。
(4.4)将CEM因截取导致的断开部分补全,重建DECNID
1)从CEM中取任一非-∞元素d1i1j生成一条边CE1=[id1i,id1j,d1i1j],其中id1i,id1j表示d1i1j对应CEM的行序号和列序号,令EXI={id1i,id1j},创建空队列QU2并将id1i,id1j入队列。
2)如QU2为空则进入步骤3),否则顶点唯一标识idi出队列QU2,在CEM中查找与idi在同行和同列的排除EXI元素的列序号和行序号idj,…idk及其非-∞元素dj,…dk,组成新建矩阵然后令/>EXI=EXI∪{idj,…idk},将idj,…idk入队列QU2,将CE1i按行拼接加入CE1生成边矩阵CE1,跳转至本步骤开头。
3)对CEM执行步骤1)至步骤2),每次执行生成新的边矩阵CE2,…CEn,直到CEM中所有元素均为-∞为止,得到边矩阵集合CE={CE1,CE2,…CEn}。
4)从距离矩阵CDM中查找出CE里面各断开部分之间的最小距离,得到边矩阵其中,ICE的第三列表示CE里面每两个边矩阵的顶点唯一标识集合(边矩阵第一二列的并集)之间的最小距离,而第一二列分别为其对应的CE中各边矩阵的序号,即边矩阵唯一标识;
5)从ICE中查找第三列的最小值,将其所在边从ICE删除后加入空的边矩阵CEn+1
6)不断从ICE中查找一个边矩阵唯一标识属于CEn+1同时另一边矩阵唯一标识属于ICE边矩阵唯一标识集合减去CEn+1边矩阵唯一标识集合的所有边,再将其中第三列最小值的边从ICE删除后按行拼接加入CEn+1,直到ICE的边矩阵唯一标识集合包含于CEn+1的边矩阵唯一标识集合为止,得到CE中各边矩阵间断开部分的补全边矩阵CEn+1
7)按行拼接边矩阵CE1,CE2,…CEn,CEn+1得到CE2=[CE1;CE2;…CEn;CEn+1],取CE2的第三列为向量DECNID
(4.5)执行步骤1中的(1.5)至(1.7)步骤。
(4.6)如CTI≤0.5则将NID(row)从NID删除再加入到CIDcol并转到步骤(4.2)继续执行,否则得到最终聚类结果簇集合CID包括密度峰值与簇边界点,以及噪点集合NID,结束所有步骤。
将上述聚类的流程步骤应用于多源大气数据聚类,建立基于分布密度的多源大气数据聚类模型如图1。图1中对于气象监测数据的定义为某一气象站在某一天时间区间内监测到的所有气象状态监测记录。气象监测数据集DS1是气象数据记录的集合,一条数据记录是某一气象站在某一时刻监测到的天气状态。一条记录包含天气、温度、气压、风速等M1个指标,其中N1为数据集中数据记录的总数。气象监测数据集DS1表示如下:
DS1={[weather,temperature,humidity,pressure,wind_direction,wind_speed]i,1≤i≤N1} (4)
图1中对于空气质量数据的定义为某一监测站在某一天时间区间内监测到的所有空气质量监测记录。空气质量数据集DS2是空气质量监测记录的集合,一条监测记录是某一监测站在某一时刻监测到的空气质量。一条记录包含PM2.5、PM10、SO2等M2个指标,其中N2为数据集中监测记录的总数。空气质量数据集DS2表示如下:
DS2={[PM2.5,PM10,NO2,CO,O3,SO2]i,1≤i≤N2} (5)
首先,合并上述数据集DS1和DS2,得到数据集DS如下(id为顶点唯一标识):
考虑到现实中存在M1≠M2,N1≠N2的情况,因此距离矩阵DM采用动态时间规整算法(Dynamic Time Warping,DTW,详见“Dynamic Programming Algorithm Optimization forSpoken Word Recognition”)计算得到DS中每两条记录的DTW距离矩阵如下:
式(7)为N×N对称矩阵,其中N=N1+N2
然后,对数据集DS和距离矩阵DM使用ANDD算法,其中涉及到顶点间距离计算的采用DTW距离,涉及到维度的采用M=M1+M2,得到包含子集合(一个子集合为一个簇)的簇集合CID与噪点集合NID,或者“无聚类趋势”提示。
最后,如果“无聚类趋势”则源1和源2的所有数据均为一天内常见天气现象出现时刻,否则按照顶点(记录)的唯一标识划分出DS1和DS2的包含记录集合,即{CID1,NID1}和{CID2,NID2}为气象监测数据和空气质量数据的聚类结果,其示意如图1,线框内为划分后的簇集合{CID1,CID2},即源1和源2的一天内常见天气现象出现时刻,线框外则为分离出的噪点集合{NID1,NID2},即源1和源2的一天内灾害性天气现象出现时刻。此外,根据该聚类结果还可以进一步的给各条记录添加类别标签作为记录的新指标来增强数据间的关联性,从而帮助后续的各种神经网络模型更好地捕捉数据间的深层特征,提高模型的灾害性天气现象预测准确率。
在本实施方式中使用三个数据集进行分析,其中数据集1为自定义数据集eye,含3050个2维数据;数据集2为flame数据集,240个2维数据;数据集3为s2数据集,5000个2维数据。
图2的eye数据集特点是具有3个大小不均的不规则簇(簇1有2500个数据、簇2和簇3有250个数据),且还有分布不均的噪点(50个噪点)。ANDD算法准确的将eye数据集的噪点分离出来并聚类了簇1至簇3,可见其具有聚类任意分布形状、分布均匀性数据的功能,并可以达到无参数自动化分离噪点的目标。
图3至图12的KNN-DPC、DPC、AP(Affinity Propagation)、DBSCAN四种算法在flame、s2数据集聚类效果的参数设置如附图说明。从图3上来看ANDD将左上角的2个噪点准确的分离了出来而其余四种算法均未能分离该噪点,其中,如图4所示,KNN-DPC未能分离噪点的原因在于其对噪点的判定是依赖于近邻样本数K参数值的,因此当调节K值改变簇的聚类结果的同时噪点的判定结果也会同步发生变化,而这却并不是希望达到的效果,此外如图5所示,DPC不具备分离噪点的功能,如图6和图7所示AP和DBSCAN则对于簇的划分结果不如KNN-DPC和DPC准确。s2数据集也是含有噪点和15个自然簇的部分簇高度重叠的数据集,从聚类效果上来看如图8所示,ANDD仍然是最好的,如图9和图10所示,KNN-DPC和DPC均未能分离出噪点且KNN-DPC还对第14个簇周围的部分离散点进行了错误分配,而如图11所示,AP只聚类出了14个簇,如图12所示,DBSCAN则正确聚类了9个簇并将离散点错误识别为第8和第13这2个簇。

Claims (5)

1.一种基于分布密度的多源大气数据聚类方法,其特征在于,包括以下步骤:
(1)构建多源大气数据集DS,并判断多源大气数据集DS的聚类趋势;所述多源大气数据集DS为一个M维数据组成的数据量为N的数据集;
(2)生成多源大气数据集的距离矩阵DM的全邻域分布密度矩阵DDM;
(3)以分布密度阈值ddth为参数,划分出全邻域分布密度矩阵DDM的密度峰值与离散点;
(4)截取全体数据的边矩阵E将簇边界点合并入簇核,得到聚类结果;
所述步骤(3)包括以下步骤:
(31)设待标记顶点唯一标识集合AID=ID,密度峰值的顶点唯一标识集合离散点的顶点唯一标识集合/>并创建空队列QU;
(32)从AID中选择第1个顶点唯一标识AID(1)并查找其在DDM中的对应行DDMrAID(1),将DDMrAID(1)按IrAID(1)重新排序得到顺序排列的DDMsrAID(1),然后从第2个元素开始依序遍历DDMsrAID(1)找到前ci个大于阈值ddth的元素,如果ci=0则将AID(1)从AID删除再加入到NID集合中并跳转至步骤(35),否则将ci个元素在DDM中对应列的唯一标识ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))和AID取交集后入队列QU,并将AID(1),ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))从AID删除再加入到新建集合CIDnew={AID(1),ID(IrAID(1)(2))…ID(IrAID(1)(ci+1))};
(33)如QU为空则进入步骤(34),否则顶点唯一标识ID(i)出队列QU,查找ID(i)在DDM中的对应行DDMrID(i),将DDMrID(i)按IrID(i)重新排序得到DDMsrID(i),然后从第2个元素开始依序遍历DDMsrID(i)找到前ci个大于阈值ddth的元素,如果ci=0跳转至本步骤开头,否则将ci个元素在DDM中对应列的唯一标识ID(IrID(i)(2))…ID(IrID(i)(ci+1))和AID取交集后入队列QU,并将ID(i),ID(IrID(i)(2))…ID(IrID(i)(ci+1))从AID删除再更新集合CIDnew=CIDnew∪{ID(i),ID(IrID(i)(2))…ID(IrID(i)(ci+1))},跳转至本步骤开头;
(34)检查最新生成并更新完成的集合CIDnew,如果CIDnew和CID集合中的某些子集合CIDj,…CIDk的交集不为空,则CIDnew=CIDnew∪CIDj∪…CIDk,并将CIDj,…CIDk从CID删除后将CIDnew加入到CID成为新增的子集合;
(35)如AID不为空则转到步骤(32)继续执行,否则得到密度峰值顶点集合CID和离散顶点集合NID并进入步骤(4)。
2.根据权利要求1所述的一种基于分布密度的多源大气数据聚类方法,其特征在于,所述步骤(1)包括以下步骤:
(11)数据集xi=[xi1...xiM]i∈{1,...N}表示数据集的第i个顶点,idi=i表示第i个顶点的唯一标识;计算距离矩阵/>其中D(xi,xj)i,j∈{1,...N}表示两个顶点间的距离,以及距离矩阵的另一种表示形式其中i,j∈{1,...N}且i≠j表示第i个和第j个顶点的唯一标识,D(xi,xj)i,j∈{1,...N}且i≠j表示两个顶点间的距离;
(12)从DE中查找第三列的最小值,将其所在行,即边,从DE删除后加入空的边矩阵E;
(13)不断从DE中查找一个顶点唯一标识属于E同时另一顶点唯一标识属于DE顶点唯一标识集合减去E顶点唯一标识集合的所有边,再将其中第三列最小值的边从DE删除后按行拼接加入E,直到DE的顶点唯一标识集合包含于E的顶点唯一标识集合为止,得到最小距离连通的边矩阵E;
(14)生成最小距离连通图MDG[ID,E],其中ID={id1,id2,…idN},然后取边矩阵E的第三列为DECNID
(15)计算DECNID中除第一行外每一行和上一行的差值,取其中的最大值Dmax在DECNID中对应行的值和对应行下一行的值计算平均值,得到均值deth,deth即为肘阈值;
(16)用DECNID中小于deth的值建立距离集合Dsmall,用DECNID中大于deth的值建立距离集合Dbig;计算出距离集合Dsmall、Dbig中元素的个数numsmall、numbig;均值meansmall、meanbig,以及Dsmall、DECNID的变异系数cvsmall、cvECNID
(17)根据以下公式得到聚类趋势指数CTI:
式中,
(18)如CTI>0.5则进入步骤(2),否则提示“无聚类趋势”信息并结束本次聚类流程。
3.根据权利要求1所述的一种基于分布密度的多源大气数据聚类方法,其特征在于,所述步骤(2)包括以下步骤:
(21)以距离矩阵DM的各元素为半径r,通过超球体公式和维度M得到超球体体积矩阵VM;
(22)对体积矩阵VM的第1行VMr1按由小到大顺序排序得到VMsr1以及VMr1元素顺序在排序后的重新排列Ir1,然后按照i={1,2…N}的顺序遍历VMsr1的体积值VMsr1(i),计算累计量同时生成VMsr1(i)对应的顶点数量PNsr1(i)=i,但如果遇到VMsr1(i)=VMsr1(i-1)=…=VMsr1(i-j)的情况则PNsr1(i)=PNsr1(i-1)=…=PNsr1(i-j)=i;
(23)将向量AVMsr1,PNsr1按照VMr1的元素顺序Ir1重新排回原来顺序得到AVMr1,PNr1
(24)对体积矩阵VM的第2至N行VMr2,VMr3…VMrN执行步骤(22)至步骤(23),得到Ir2,…IrN和AVMr2,…AVMrN和PNr2,…PNrN
(25)将向量集合{Ir1,Ir2,…IrN}依序逐行组合成行元素顺序排序矩阵I,I为N×N非对称阵,将向量集合{AVMr1,AVMr2,…AVMrN}依序逐行组合成超球体累计体积矩阵AVM,AVM为N×N非对称阵,将向量集合{PNr1,PNr2,…PNrN}依序逐行组合成顶点数量矩阵PN,PN为N×N非对称阵,计算得到DDM,DDM为N×N非对称矩阵如下:
4.根据权利要求1所述的一种基于分布密度的多源大气数据聚类方法,其特征在于,所述步骤(4)包括以下步骤:
(41)创建N×N矩阵EM赋值-∞,然后遍历N个顶点的MDG[ID,E]中边矩阵E的各行,即各条边eqp=[idq,idp,dqp]令得到E的另一种表示形式EM;
(42)从距离矩阵DM中截取离散顶点集合NID为行序号,密度峰值顶点集合CID子集合的并集为列序号的部分,求其最小值,得到其中NID(row)表示最小值对应的离散顶点唯一标识,CID(col)表示最小值对应的密度峰值顶点唯一标识,而其所属的CID子集合表示为CIDcol
(43)从DM、EM中截取CNID=CIDcol∪NID(row)部分,得到CDM=DM(CNID,CNID)、CEM=EM(CNID,CNID);取CEM中非-∞元素组成向量DECNID,如果DECNID的元素数量小于CNID的元素数量减一,则进入步骤(44),否则进入步骤(45);
(44)将CEM因截取导致的断开部分补全,重建DECNID
(45)执行步骤(15)至(17);
(46)如CTI≤0.5则将NID(row)从NID删除再加入到CIDcol并转到步骤(42)继续执行,否则得到最终聚类结果簇集合CID包括密度峰值与簇边界点,以及噪点集合NID,结束所有步骤。
5.根据权利要求4所述的一种基于分布密度的多源大气数据聚类方法,其特征在于,所述步骤(44)包括以下步骤:
(441)从CEM中取任一非-∞元素d1i1j生成一条边CE1=[id1i,id1j,d1i1j],其中id1i,id1j表示d1i1j对应CEM的行序号和列序号,令创建空队列QU2并将id1i,id1j入队列;
(442)如QU2为空则进入步骤(443),否则顶点唯一标识idi出队列QU2,在CEM中查找与idi在同行和同列的排除EXI元素的列序号和行序号idj,…idk及其非-∞元素dj,…dk,组成新建矩阵然后令/>EXI=EXI∪{idj,…idk},将idj,…idk入队列QU2,将CE1i按行拼接加入CE1生成边矩阵CE1,跳转至本步骤开头;
(443)对CEM执行步骤(441)至步骤(442),每次执行生成新的边矩阵CE2,…CEn,直到CEM中所有元素均为-∞为止,得到边矩阵集合CE={CE1,CE2,…CEn};
(444)从距离矩阵CDM中查找出CE里面各断开部分之间的最小距离,得到边矩阵其中ICE的第三列表示CE里面每两个边矩阵的顶点唯一标识集合(边矩阵第一二列的并集)之间的最小距离,而第一二列分别为其对应的CE中各边矩阵的序号,即边矩阵唯一标识;
(445)从ICE中查找第三列的最小值,将其所在边从ICE删除后加入空的边矩阵CEn+1
(446)不断从ICE中查找一个边矩阵唯一标识属于CEn+1同时另一边矩阵唯一标识属于ICE边矩阵唯一标识集合减去CEn+1边矩阵唯一标识集合的所有边,再将其中第三列最小值的边从ICE删除后按行拼接加入CEn+1,直到ICE的边矩阵唯一标识集合包含于CEn+1的边矩阵唯一标识集合为止,得到CE中各边矩阵间断开部分的补全边矩阵CEn+1
(447)按行拼接边矩阵CE1,CE2,…CEn,CEn+1得到CE2=[CE1;CE2;…CEn;CEn+1],取CE2的第三列为向量DECNID
CN202010314605.3A 2020-04-21 2020-04-21 一种基于分布密度的多源大气数据聚类方法 Active CN111507415B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010314605.3A CN111507415B (zh) 2020-04-21 2020-04-21 一种基于分布密度的多源大气数据聚类方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010314605.3A CN111507415B (zh) 2020-04-21 2020-04-21 一种基于分布密度的多源大气数据聚类方法

Publications (2)

Publication Number Publication Date
CN111507415A CN111507415A (zh) 2020-08-07
CN111507415B true CN111507415B (zh) 2023-07-25

Family

ID=71876289

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010314605.3A Active CN111507415B (zh) 2020-04-21 2020-04-21 一种基于分布密度的多源大气数据聚类方法

Country Status (1)

Country Link
CN (1) CN111507415B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112232590B (zh) * 2020-11-02 2023-06-30 国网湖南省电力有限公司 多源电力气象融合数据整体评估方法、系统及存储介质
CN112904276B (zh) * 2021-01-25 2022-11-25 中国气象科学研究院 闪电辐射源连接方法
CN112734777B (zh) * 2021-01-26 2022-10-11 中国人民解放军国防科技大学 一种基于簇形状边界闭包聚类的图像分割方法及系统
CN113158817B (zh) * 2021-03-29 2023-07-18 南京信息工程大学 一种基于快速密度峰聚类的客观天气分型方法
CN112990355A (zh) * 2021-04-15 2021-06-18 中科三清科技有限公司 污染天气的分型方法、装置、电子设备及存储介质
CN117743876A (zh) * 2023-12-22 2024-03-22 冻冻(北京)网络科技有限公司 基于云计算的智慧仓储数据优化管理方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107340364A (zh) * 2017-05-31 2017-11-10 北京市环境保护监测中心 基于海量大气污染浓度数据的污染空间分析方法及装置
CN108448610A (zh) * 2018-03-12 2018-08-24 华南理工大学 一种基于深度学习的短期风功率预测方法
CN109978228A (zh) * 2019-01-31 2019-07-05 中南大学 一种pm2.5浓度预测方法、装置及介质
CN110245692A (zh) * 2019-05-27 2019-09-17 南京信息工程大学 一种用于集合数值天气预报成员的层次聚类方法
CN110263814A (zh) * 2019-05-27 2019-09-20 南京信息工程大学 基于动态聚类趋势分析的增量聚类数据挖掘方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107340364A (zh) * 2017-05-31 2017-11-10 北京市环境保护监测中心 基于海量大气污染浓度数据的污染空间分析方法及装置
CN108448610A (zh) * 2018-03-12 2018-08-24 华南理工大学 一种基于深度学习的短期风功率预测方法
CN109978228A (zh) * 2019-01-31 2019-07-05 中南大学 一种pm2.5浓度预测方法、装置及介质
CN110245692A (zh) * 2019-05-27 2019-09-17 南京信息工程大学 一种用于集合数值天气预报成员的层次聚类方法
CN110263814A (zh) * 2019-05-27 2019-09-20 南京信息工程大学 基于动态聚类趋势分析的增量聚类数据挖掘方法

Also Published As

Publication number Publication date
CN111507415A (zh) 2020-08-07

Similar Documents

Publication Publication Date Title
CN111507415B (zh) 一种基于分布密度的多源大气数据聚类方法
CN107808518B (zh) 一种基于多路谱聚类理论的交通小区划分方法
CN111832655B (zh) 一种基于特征金字塔网络的多尺度三维目标检测方法
CN109753995B (zh) 一种基于PointNet++的3D点云目标分类和语义分割网络的优化方法
WO2012111236A1 (ja) 画像識別装置およびプログラム
CN108614997B (zh) 一种基于改进AlexNet的遥感图像识别方法
WO2016066042A1 (zh) 商品图片的分割方法及其装置
CN105354593B (zh) 一种基于nmf的三维模型分类方法
CN104166982A (zh) 基于典型相关性分析的图像优化聚类方法
CN106203494B (zh) 一种基于内存计算的并行化聚类方法
CN110751027B (zh) 一种基于深度多示例学习的行人重识别方法
CN106780639B (zh) 基于显著性特征稀疏嵌入和极限学习机的哈希编码方法
CN110866896A (zh) 基于k-means与水平集超像素分割的图像显著性目标检测方法
CN106845536B (zh) 一种基于图像缩放的并行聚类方法
CN101339553A (zh) 面向海量数据近似快速聚类和索引方法
CN110047139A (zh) 一种指定目标三维重建方法及系统
CN110543895A (zh) 一种基于VGGNet和ResNet的图像分类方法
CN110349159A (zh) 基于权重能量自适应分布的三维形状分割方法及系统
CN103617609A (zh) 基于图论的k-means非线性流形聚类与代表点选取方法
CN104732524A (zh) 血液白细胞显微图像的随机权网络分割方法
CN113222181A (zh) 一种面向k-means聚类算法的联邦学习方法
CN110084136A (zh) 基于超像素crf模型的上下文优化室内场景语义标注方法
CN116993555A (zh) 国土空间规划重点区域识别的分区方法、系统及存储介质
CN112711985B (zh) 基于改进solo网络的果实识别方法、装置及机器人
Kumar et al. A hybrid cluster technique for improving the efficiency of colour image segmentation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant