CN106127809B - 一种显微图像序列中癌细胞轨迹追踪与关联方法 - Google Patents

一种显微图像序列中癌细胞轨迹追踪与关联方法 Download PDF

Info

Publication number
CN106127809B
CN106127809B CN201610458726.9A CN201610458726A CN106127809B CN 106127809 B CN106127809 B CN 106127809B CN 201610458726 A CN201610458726 A CN 201610458726A CN 106127809 B CN106127809 B CN 106127809B
Authority
CN
China
Prior art keywords
cell
track
frame
cancer cell
state
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
CN201610458726.9A
Other languages
English (en)
Other versions
CN106127809A (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.)
Zhejiang University of Technology ZJUT
Original Assignee
Zhejiang University of Technology ZJUT
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 Zhejiang University of Technology ZJUT filed Critical Zhejiang University of Technology ZJUT
Priority to CN201610458726.9A priority Critical patent/CN106127809B/zh
Publication of CN106127809A publication Critical patent/CN106127809A/zh
Application granted granted Critical
Publication of CN106127809B publication Critical patent/CN106127809B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10016Video; Image sequence
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10056Microscopic image
    • G06T2207/10061Microscopic image from scanning electron microscope
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing

Landscapes

  • Investigating Or Analysing Biological Materials (AREA)
  • Image Analysis (AREA)

Abstract

一种显微图像序列中癌细胞细胞轨迹追踪与关联方法,包括如下步骤:1)基于最近邻法的癌细胞轨迹片段形成,将检测到的癌细胞区域前后帧关联,形成可信的轨迹片段;2)基于细胞动态特征的全局轨迹关联;3)基于分类检测的轨迹关联;4)基于轨迹关联的状态检测纠正,在癌细胞生命轨迹关联之后,即观测细胞都已形成连续的轨迹后,根据前后帧细胞的状态对于误检测的细胞状态进行修正。本发明既能实现对癌细胞细胞轨迹的关联,同时利用分类检测和轨迹关联的状态检测又能更好的关联轨迹片段及对其加以修正。

Description

一种显微图像序列中癌细胞轨迹追踪与关联方法
技术领域
本发明涉及细胞成像技术研究领域,尤其是一种癌细胞细胞轨迹跟踪与关联方法。
背景技术
细胞轨迹追踪就是结合视频前后帧的信息以及细胞分裂生长的过程,实现对连续帧间细胞生命轨迹的追踪和关联。而本发明的目标对象是受药物抑制的癌细胞,其分裂生长过程较正常细胞具有特殊性,复杂性。传统的单一的细胞追踪方法并不适用。
从目前已有的多种目标跟踪方法来看,近年来研究较多的多目标的细胞跟踪主要分为两类。第一类基于帧间数据关联的算法,即通过前一帧信息来估计当前帧细胞情况,如此循环迭代。第二类是通过全局数据关联的多目标跟踪方法,利用更大的时间窗口获得更多的信息,通过全局轨迹长期的检测,在给定时间段内全局的解决轨迹连接问题。这些方法在处理树形层级轨迹关联时(即母细胞分裂到子细胞的情况)仍有一定的问题,需要我们在细胞分类后连续观测几帧细胞的情况以对细胞关联关系进行修正。
发明内容
为了解决现有技术中一个目标对象的跟踪门限中量测数据可能不止一个(即量测数据可能来自于正确目标对象,也有可能来自于其他目标对象,也可能来自于杂波)的问题、精度较低,本发明提出了一种精度较高的显微图像序列中癌细胞轨迹追踪与关联方法。
本发明解决其技术问题所采用的技术方案是:
一种显微图像序列中癌细胞轨迹追踪与关联方法,所述细胞轨迹追踪与关联方法包括如下步骤:
1)基于最近邻法的癌细胞轨迹片段形成,即将检测到的癌细胞区域前后帧关联,形成可信的轨迹片段,过程如下:
首先,设置跟踪波门限r,落到跟踪波门中的量测作为候选回波,即目标的量测值zi(k)是否能够满足公式:
其中是跟踪波门的中心,如果波门中只有一个量测,该量测值则直接被用于轨迹的更新中,形成连续的轨迹段;如若有一个以上的候选回波,则应根据距离计算方式,找到距离最近的候选回波,来用于轨迹的更新;其中,zi(k)对应的统计距离为:
接着,考虑当跟踪波门中只有一个量测,或是距离最近的候选回波距离前一帧目标的距离d足够近的情况下,将其连接到细胞轨迹中,形成细胞轨迹片段;而当最近量测值距离超过某一阈值,或是有两个及以上候选回波距离前一帧目标的距离d都比较近的情况下,细胞轨迹片段断裂,重新开始新的轨迹;
2)基于细胞动态特征的全局轨迹关联;
由于细胞在生长、分裂、融合的过程中形态不断发生改变,仅依靠单一的特征难以正确区分,因此结合细胞的动态特征匹配,将再次断裂的细胞轨迹片段关联。
细胞动态特征参数包括以下4种:运动位移参数Edisplacement、运动偏斜参数Eskewness、面积变化参数Earea以及形状变化参数Edeformation
细胞u和v的运动位移参数Edisplacement定义如下:
其中,cu,cv分别是细胞u、v在第n帧和第n+1帧的位置,H和W分别指当前帧图像的长和宽;
运动偏斜参数Eskewness定义如下:
这里ct,cu,cv分别是细胞t,u,v在第n-1帧,第n帧以及第n+1帧时的中心位置;这个参数是用了测量运动方向的偏移。
面积变化参数Earea用于测量细胞u和v的重叠程度,定义如下:
其中,指细胞u在第n帧的面积,指细胞v在第n+1帧中的面积,指两细胞的重叠程度;这个公式表明,两细胞重叠面积越多,两者就更倾向于是同一个细胞。
对于形状变化参数参数Edeformation,使用匹配细胞u和v面积的椭圆离心率,椭圆离心率公式Q=P2/(4π×A2),这里P和A分别表示椭圆的周长和面积,形变参数Edeformation定义如下:
使用4个参数来估计轨迹片段间的相似程度,另外,使用一个基于双边权值匹配的全局最优技术来最小化总代价,函数定义如下:
E(ci,cj)=α1Edisplacement2Eskewness3Earea4Edeformation
其中,E(ci,cj)就是第n帧中的细胞i与第n+1帧中的细胞j的代价函数,αi限定在0到1之间,且
通过计算各个细胞的动态特征,得到代价函数的值之后,值最小的就是正确匹配的细胞,即可将此处断裂的细胞轨迹片段正确关联起来。
3)基于分类检测的轨迹关联,过程如下:
在数据关联中,无论是基于最近邻法的癌细胞轨迹片段形成还是基于动态特征的全局轨迹关联,都是一对一的匹配过程。如若涉及到细胞分裂的情况,则需要根据分裂的具体状况进行分类检测以关联轨迹。
首先,根据细胞分裂时候的特征判断是否属于分裂事件;
再者,因为本发明研究对象为受药物抑制的癌细胞,正常细胞的分裂检测事件判断对本数据集并不适用。根据癌细胞分裂特征,判断是否将未关联的轨迹片段关联到母细胞,并标记分裂事件。
4)基于轨迹关联的状态检测纠正;
在癌细胞生命轨迹关联之后,即观测细胞都已形成连续的轨迹后,根据前后帧细胞的状态对于误检测的细胞状态进行修正;
从第一帧开始,记录下待监测癌细胞所处的阶段,在该细胞的整个轨迹过程中,向后逐步检测,记录下与前一帧状态不同的时刻;若所记录的时刻为突变点,即相对整个序列的异常点,则根据前后细胞状态修正误检测的情况;
另外,由于细胞生长遵循分裂间期、分裂中期和分裂后期这三个阶段,细胞轨迹中的状态一定也是按此顺序出现,由于各个癌细胞所处阶段不同,起始状态可能是分裂间期、分裂前期或是分裂中期,但在分裂后期之前,都遵循正常的规律。分裂后期之后,即受药物抑制癌细胞挣扎后的状态可能会有所不同;
遵循如上两个原则以修正癌细胞细胞跟踪轨迹。
本发明的有益效果主要表现在:基于数据关联的癌细胞细胞轨迹追踪方法,既能实现对癌细胞细胞轨迹的关联,同时利用分类检测和轨迹关联的状态检测又能更好的关联轨迹片段及对其加以修正。
具体实施方式
下面对本发明作进一步描述。
一种显微图像序列中癌细胞轨迹追踪与关联方法,所述追踪与关联方法包括如下步骤:
1)基于最近邻法的癌细胞轨迹片段形成,即将检测到的癌细胞区域前后帧关联,形成可信的轨迹片段,过程如下:
首先,设置跟踪波门限r,落到跟踪波门中的量测作为候选回波,即目标的量测值zi(k)是否能够满足公式:
其中是跟踪波门的中心,如果波门中只有一个量测,该量测值则直接被用于轨迹的更新中,形成连续的轨迹段;如若有一个以上的候选回波,则应根据距离计算方式,找到距离最近的候选回波,来用于轨迹的更新;其中,zi(k)对应的统计距离为:
接着,考虑当跟踪波门中只有一个量测,或是距离最近的候选回波距离前一帧目标的距离d足够近的情况下,将其连接到细胞轨迹中,形成细胞轨迹片段。而当最近量测值距离超过某一阈值,或是有两个及以上候选回波距离前一帧目标的距离d都比较近的情况下,细胞轨迹片段断裂,重新开始新的轨迹;
2)基于细胞动态特征的全局轨迹关联;
由于细胞在生长、分裂、融合的过程中形态不断发生改变,仅依靠单一的特征难以正确区分,因此结合细胞的动态特征匹配,将再次断裂的细胞轨迹片段关联;
细胞动态特征参数包括以下4种:运动位移参数(Edisplacement)、运动偏斜参数(Eskewness)、面积变化参数(Earea)、以及形状变化参数(Edeformation)。
细胞u和v的运动位移参数Edisplacement定义如下:
其中,cu,cv分别是细胞u、v在第n帧和第n+1帧的位置,H和W分别指当前帧图像的长和宽。
运动偏斜参数Eskewness定义如下:
其中,ct,cu,cv分别是细胞t,u,v在第n-1帧,第n帧以及第n+1帧时的中心位置,运动偏斜参数是用了测量运动方向的偏移;
面积变化参数Earea用于测量细胞u和v的重叠程度,定义如下:
其中,指细胞u在第n帧的面积,指细胞v在第n+1帧中的面积,指两细胞的重叠程度。这个公式表明,两细胞重叠面积越多,两者就更倾向于是同一个细胞。
对于形状变化参数Edeformation,使用匹配细胞u和v面积的椭圆离心率,椭圆离心率公式Q=P2/(4π×A2),这里P和A分别表示椭圆的周长和面积,形状变化参数Edeformation定义如下:
使用4个参数来估计轨迹片段间的相似程度,另外,使用一个基于双边权值匹配的全局最优技术来最小化总代价,函数定义如下:
E(ci,cj)=α1Edisplacement2Eskewness3Earea4Edeformation
这里的E(ci,cj)就是第n帧中的细胞i与第n+1帧中的细胞j的代价函数。这里αi限定在0到1之间,且
通过计算各个细胞的动态特征,得到代价函数的值之后,值最小的就是正确匹配的细胞,即可将此处断裂的细胞轨迹片段正确关联起来;
3)基于分类检测的轨迹关联
在数据关联中,无论是基于最近邻法的癌细胞轨迹片段形成还是基于动态特征的全局轨迹关联,都是一对一的匹配过程。如若涉及到细胞分裂的情况,则需要根据分裂的具体状况进行分类检测以关联轨迹。
首先,根据细胞分裂时候的特征判断是否属于分裂事件,特征如下:
正常细胞发生分裂时,在下一刻会有两个独立的细胞产生,这两个细胞距离前一时刻的母细胞距离较近;
母细胞在一份为二之前,细胞形状变圆,周围发亮;
两个子细胞之间相距很近;
独立的两个子细胞面积相对较小,二者面积和近似于前一帧中母细胞的面积。
再者,因为本发明研究对象为受药物抑制的癌细胞,正常细胞的分裂检测事件判断对本数据集并不适用。首先,癌细胞若发生分裂,未必是一分为二,可能分裂为两个或是两个以上子细胞;其次,母细胞形状变圆,变圆发亮这个阶段,药物抑制下的癌细胞也会经历(即第二个阶段),但通常在这一阶段后癌细胞会与药物对抗做挣扎,若有分裂也是在挣扎失败之后发生,因此变圆发亮这一特征无法作为判别依据;此外,癌细胞在受药物作用挣扎的过程中运动剧烈,与之前相比可能有较大的位移偏差。根据癌细胞的这些分裂特征,对于在时空序列中途出现的未能与前帧关联起来的癌细胞轨迹片段在前一帧中该细胞附近寻找是否有处于第三阶段的细胞存在,若有,将未关联的轨迹片段关联到母细胞,并标记分裂事件。若未找到合适的匹配,则认为该细胞为由外部进入的新出现细胞;
4)基于轨迹关联的状态检测纠正
在癌细胞生命轨迹关联之后,即观测细胞都已形成连续的轨迹后,根据前后帧细胞的状态对于误检测的细胞状态进行修正。
从第一帧开始,记录下待监测癌细胞所处的阶段,在该细胞的整个轨迹过程中,向后逐步检测,记录下与前一帧状态不同的时刻。若所记录的时刻为突变点,即相对整个序列的异常点,则根据前后细胞状态修正误检测的情况。
另外,由于细胞生长遵循分裂间期、分裂中期和分裂后期这三个阶段,细胞轨迹中的状态一定也是按此顺序出现。由于各个癌细胞所处阶段不同,起始状态可能是分裂间期、分裂前期或是分裂中期,但在分裂后期之前,都遵循正常的规律。分裂后期之后,即受药物抑制癌细胞挣扎后的状态可能会有所不同。
遵循如上两个原则,选取异常检测方法,如常见的基于统计的方法,基于距离的方法,基于密度的方法,基于偏离的方法,基于聚类的方法等以修正癌细胞细胞跟踪轨迹。

Claims (3)

1.一种显微图像序列中癌细胞轨迹追踪与关联方法,其特征在于:所述细胞轨迹追踪与关联方法包括如下步骤:
1)基于最近邻法的癌细胞轨迹片段形成,即将检测到的癌细胞区域前后帧关联,形成可信的轨迹片段,过程如下:
首先,设置跟踪波门限r,落到跟踪波门中的量测作为候选回波,即目标的量测值zi(k)是否能够满足公式:
其中,是跟踪波门的中心,如果波门中只有一个量测,该量测值则直接被用于轨迹的更新中,形成连续的轨迹段;如若有一个以上的候选回波,则应根据距离计算方式,找到距离最近的候选回波,来用于轨迹的更新;其中,zi(k)对应的统计距离为:
接着,考虑当跟踪波门中只有一个量测,或是距离最近的候选回波距离前一帧目标的距离d足够近的情况下,将其连接到细胞轨迹中,形成细胞轨迹片段;而当最近量测值距离超过某一阈值,或是有两个及以上候选回波距离前一帧目标的距离d都比较近的情况下,细胞轨迹片段断裂,重新开始新的轨迹;
2)基于细胞动态特征的全局轨迹关联;
细胞动态特征参数包括以下4种:运动位移参数Edisplacement、运动偏斜参数Eskewness、面积变化参数Earea以及形状变化参数Edeformation
细胞u和v的运动位移参数Edisplacement定义如下:
其中,cu,cv分别是细胞u、v在第n帧和第n+1帧的位置,H和W分别指当前帧图像的长和宽;
运动偏斜参数Eskewness定义如下:
这里ct,cu,cv分别是细胞t,u,v在第n-1帧,第n帧以及第n+1帧时的中心位置;
面积变化参数Earea用于测量细胞u和v的重叠程度,定义如下:
其中,指细胞u在第n帧的面积,指细胞v在第n+1帧中的面积,指两细胞的重叠程度;
对于形状变化参数Edeformation,使用匹配细胞u和v面积的椭圆离心率,椭圆离心率公式Q=P2/(4π×A2),这里P和A分别表示椭圆的周长和面积,形变参数Edeformation定义如下:
使用一个基于双边权值匹配的全局最优技术来最小化总代价,函数定义如下:
E(ci,cj)=α1Edisplacement2Eskewness3Earea4Edeformation
其中,E(ci,cj)就是第n帧中的细胞i与第n+1帧中的细胞j的代价函数,αi限定在0到1之间,且
通过计算各个细胞的动态特征,得到代价函数的值之后,值最小的就是正确匹配的细胞,即可将此处断裂的细胞轨迹片段正确关联起来;
3)基于分类检测的轨迹关联,过程如下:
首先,根据细胞分裂时候的特征判断是否属于分裂事件;
再者,根据癌细胞分裂特征,判断是否将未关联的轨迹片段关联到母细胞,并标记分裂事件;
4)基于轨迹关联的状态检测纠正;
在癌细胞生命轨迹关联之后,即观测细胞都已形成连续的轨迹后,根据前后帧细胞的状态对于误检测的细胞状态进行修正。
2.如权利要求1所述的一种显微图像序列中癌细胞轨迹追踪与关联方法,其特征在于:所述步骤4)中,从第一帧开始,记录下待监测癌细胞所处的阶段,在该细胞的整个轨迹过程中,向后逐步检测,记录下与前一帧状态不同的时刻;若所记录的时刻为突变点,即相对整个序列的异常点,则根据前后细胞状态修正误检测的情况。
3.如权利要求2所述的一种显微图像序列中癌细胞轨迹追踪与关联方法,其特征在于:由于细胞生长遵循分裂间期、分裂中期和分裂后期这三个阶段,细胞轨迹中的状态也是按此顺序出现;由于各个癌细胞所处阶段不同,起始状态可能是分裂间期、分裂前期或是分裂中期,但在分裂后期之前,都遵循正常的规律,分裂后期之后,即受药物抑制癌细胞挣扎后的状态可能会有所不同。
CN201610458726.9A 2016-06-22 2016-06-22 一种显微图像序列中癌细胞轨迹追踪与关联方法 Active CN106127809B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610458726.9A CN106127809B (zh) 2016-06-22 2016-06-22 一种显微图像序列中癌细胞轨迹追踪与关联方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610458726.9A CN106127809B (zh) 2016-06-22 2016-06-22 一种显微图像序列中癌细胞轨迹追踪与关联方法

Publications (2)

Publication Number Publication Date
CN106127809A CN106127809A (zh) 2016-11-16
CN106127809B true CN106127809B (zh) 2018-09-18

Family

ID=57267837

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610458726.9A Active CN106127809B (zh) 2016-06-22 2016-06-22 一种显微图像序列中癌细胞轨迹追踪与关联方法

Country Status (1)

Country Link
CN (1) CN106127809B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109118492B (zh) * 2016-12-22 2022-04-05 浙江工业大学 一种相差显微镜细胞图像的检测和分割方法
CN108665483B (zh) * 2018-05-04 2021-10-12 浙江工业大学 一种基于多特征融合的癌细胞跟踪方法
CN109523577B (zh) * 2018-10-29 2020-09-01 浙江大学 基于显微图像的亚细胞结构运动轨迹确定方法
WO2023019559A1 (zh) * 2021-08-20 2023-02-23 深圳先进技术研究院 一种自动化干细胞检测方法、系统、终端以及存储介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2725549A1 (en) * 2012-10-26 2014-04-30 Varian Medical Systems Technologies, Inc. Template matching method for image-based detection and tracking of irregular shaped targets
CN103927763A (zh) * 2014-03-24 2014-07-16 河海大学 一种图像序列多目标跟踪轨迹识别处理方法
CN104915970A (zh) * 2015-06-12 2015-09-16 南京邮电大学 一种基于轨迹关联的多目标跟踪方法
CN105678810A (zh) * 2016-01-12 2016-06-15 上海交通大学 基于全局和局部层次最优的细胞跟踪方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2725549A1 (en) * 2012-10-26 2014-04-30 Varian Medical Systems Technologies, Inc. Template matching method for image-based detection and tracking of irregular shaped targets
CN103927763A (zh) * 2014-03-24 2014-07-16 河海大学 一种图像序列多目标跟踪轨迹识别处理方法
CN104915970A (zh) * 2015-06-12 2015-09-16 南京邮电大学 一种基于轨迹关联的多目标跟踪方法
CN105678810A (zh) * 2016-01-12 2016-06-15 上海交通大学 基于全局和局部层次最优的细胞跟踪方法

Also Published As

Publication number Publication date
CN106127809A (zh) 2016-11-16

Similar Documents

Publication Publication Date Title
CN106127809B (zh) 一种显微图像序列中癌细胞轨迹追踪与关联方法
Frossard et al. End-to-end learning of multi-sensor 3d tracking by detection
CN108921873B (zh) 基于核相关滤波优化的马尔科夫决策在线多目标跟踪方法
CN103942536B (zh) 一种迭代更新轨迹模型的多目标跟踪方法
CN106780557A (zh) 一种基于光流法和关键点特征的运动目标跟踪方法
CN110013247A (zh) 一种心电图p波和t波的检测、区分与定位方法
KR101402206B1 (ko) 운동학 및 형상 정보를 활용한 다중 표적 추적 방법
Houenou et al. A track-to-track association method for automotive perception systems
CN104019816A (zh) 基于全局时空相关的概率假设密度滤波器航迹提取方法
CN109118514B (zh) 一种目标跟踪方法
CN111781608A (zh) 一种基于fmcw激光雷达的运动目标检测方法及系统
CN106291498B (zh) 一种基于粒子滤波的检测跟踪联合优化方法
Damen et al. Recognizing linked events: Searching the space of feasible explanations
CN113537077B (zh) 基于特征池优化的标签多伯努利视频多目标跟踪方法
CN113671481A (zh) 一种基于毫米波雷达的3d多目标跟踪处理方法
CN103646242A (zh) 基于最大稳定极值区域特征的扩展目标跟踪方法
CN102034246A (zh) 基于拓扑约束和匈牙利算法的高密度细胞追踪方法
CN111308459A (zh) 一种基于传感器特性和量测来源的多传感器航迹起始方法
CN110363791B (zh) 一种融合单目标跟踪结果的在线多目标跟踪方法
CN113012193A (zh) 一种基于深度学习的多行人跟踪方法
CN108828584B (zh) 基于航迹折叠因子解模糊的多重频目标检测前跟踪方法
JP4488185B2 (ja) 目標追尾装置
CN107561506A (zh) 一种多级虚警概率联合递推目标检测方法
Su et al. Real-time multiple object tracking based on optical flow
Butler et al. In defense of Kalman filtering for polyp tracking from colonoscopy videos

Legal Events

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