CN103489193B - 基于融合策略的面向对象的高分辨率遥感影像变化检测方法 - Google Patents

基于融合策略的面向对象的高分辨率遥感影像变化检测方法 Download PDF

Info

Publication number
CN103489193B
CN103489193B CN201310460801.1A CN201310460801A CN103489193B CN 103489193 B CN103489193 B CN 103489193B CN 201310460801 A CN201310460801 A CN 201310460801A CN 103489193 B CN103489193 B CN 103489193B
Authority
CN
China
Prior art keywords
image
segmentation
variance
resolution
shape
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.)
Expired - Fee Related
Application number
CN201310460801.1A
Other languages
English (en)
Other versions
CN103489193A (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.)
Hohai University HHU
Changjiang River Scientific Research Institute Changjiang Water Resources Commission
Original Assignee
Hohai University HHU
Changjiang River Scientific Research Institute Changjiang Water Resources Commission
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 Hohai University HHU, Changjiang River Scientific Research Institute Changjiang Water Resources Commission filed Critical Hohai University HHU
Priority to CN201310460801.1A priority Critical patent/CN103489193B/zh
Publication of CN103489193A publication Critical patent/CN103489193A/zh
Application granted granted Critical
Publication of CN103489193B publication Critical patent/CN103489193B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Image Analysis (AREA)

Abstract

本发明公开一种基于融合策略的面向对象的高分辨率遥感影像变化检测方法,输入两时相已配准的高分辨率光学遥感影像,对输入的高分辨率影像进行多分辨率分割;对分割的单元(对象)进行变化的检测,计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域;将测结果进行基于投票法的融合,以获得更加稳健和精度较高的变化检测结果。

Description

基于融合策略的面向对象的高分辨率遥感影像变化检测方法
技术领域
本发明涉及高空间分辨率光学遥感影像变化检测,具体是一种基于融合策略的面向对象的高分辨率遥感影像变化检测方法,该方法可以用来解决高空间分辨率多光谱遥感影像变化检测精度不高的问题,属于遥感影像处理技术领域。
背景技术
随着多时相高分辨率遥感数据的不断积累以及空间数据库的相继建立,如何从这些遥感数据中提取和检测变化信息已成为遥感科学和地理信息科学的重要研究课题。根据同一区域不同时相的遥感影像,可以提取城市、环境等动态变化的信息,为资源管理与规划、环境保护等部门提供科学决策的依据。我国“十二五”将加大拓展实施“十一五”已启动实施的高分辨率对地观测工程,关注包括高分辨率遥感目标与空间环境特征分析及高可靠性自动解译等基础理论与关键技术研究,正在成为解决国家安全和社会经济发展重大需求的研究焦点。
遥感影像的变化检测就是从不同时期的遥感数据中,定量地分析和确定地表变化的特征与过程。各国学者从不同的角度和应用研究提出了许多有效的检测算法,如变化矢量分析法(ChangeVectorAnalysis,CVA)、分类后比较法等。然而,上述这些传统的变化检测方法应用于高分辨率(光学)遥感影像时,又会出现新的难点和问题,比如(1)空间分辨率的提高会产生“同物异谱”和“同谱异物”现象,从而导致变化类和非变化类的可分性降低;(2)图像相邻像素之间高度相关。
针对上述难点和存在的问题,许多学者试图通过对传统算法的改进,或者引入诸(较)多约束条件。尽管这些检测方法都具有一定的潜力,但是高分辨率影像检测环境的复杂化以及目标先验信息匮乏等,导致这些算法都存在着一定的局限性。为此,有必要研究新的高分辨率可见光遥感图像变化检测技术来有效克服上述难点。
发明内容
发明目的:高空间分辨率遥感影像背景信息复杂、噪声干扰严重,使得变化检测所面临的问题采用常规的变化检测方法难以解决。本发明在对高分辨率遥感影像进行分割的基础上,对分割的单元(对象)进行变化的检测,并采用融合的策略,将多种检测结果进行基于投票法的融合,以获得更加稳健和精度较高的变化检测结果。
技术方案:一种基于融合策略的面向对象的高分辨率遥感影像变化检测方法,包括下述步骤:
步骤1:输入两时相已配准的高分辨率光学遥感影像,分别记为:X1和X2
步骤2:对输入的高分辨率影像进行多分辨率分割:
2a)分别在尺度l1、在尺度l2、在尺度l3上分别对X1、X2、X1和X2组合的影像进行多尺度分割,分割后的影像分别记为S1、S2和S3
2b)将S1叠加到X2,根据S1的每一个对象边界确定X2的相对应的对象边界,从而获得X2的分割影像X2S1
2c)将S2叠加到X1,根据S2的每一个对象边界确定X1的相对应的对象边界,从而获得X1的分割影像X1S2
2d)将S3分别叠加到X1和X2,根据S3的每一个对象边界分别确定X1和X2的相对应的对象边界,从而获得X1和X2的分割影像X1S3和X2S3
2e)将S1和S2进行并集操作,获得分割影像S12
2f)将S12分别叠加到X1和X2,根据S12的每一个对象边界分别确定X1和X2的相对应的对象边界,从而获得X1和X2的分割影像X1S12和X2S12
步骤3:对所得到的S1和X2S1进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域。
步骤4:对所得到的S2和X1S2进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域。
步骤5:对所得到的X1S3和X2S3进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域。
步骤6:对所得到的X1S12和X2S12进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域。
步骤7:将步骤3-6的检测结果进行基于投票法的融合,获得最终的检测结果。
有益效果:与现有技术相比,本发明提供的基于融合策略的面向对象的高分辨率遥感影像变化检测方法,具有如下优点:
(1)在变化检测中采用面向对象的方法,可以提取更多对象特征,特别是对象间的空间关系,使得遥感的变化检测更有利于知识的结合应用,也更能有效地利用影像所具有的多特征优势。
(2)变化检测中,采用融合的策略,融合多种方法的检测结果,可以使得变化检测的结果更加可靠,也更加具有稳健性。
附图说明
图1为本发明方法实现流程图。
具体实施方式
下面结合具体实施例,进一步阐明本发明,应理解这些实施例仅用于说明本发明而不用于限制本发明的范围,在阅读了本发明之后,本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
如图1所示,基于融合策略的面向对象的高分辨率遥感影像变化检测方法,实现步骤如下:
步骤1:输入两时相已配准的高分辨率光学遥感影像,分别记为:X1和X2
步骤2:对输入的高分辨率影像进行多分辨率分割,在尺度hsc下,采用的多分辨率分割包括如下步骤:
1)计算对象的光谱(spectral)异质性指标hcolor
h color = Σ c w c · σ c - - - ( 1 )
其中ωc为图层的权重,σc为图层的标准差,c为图层数;根据不同的影像特性以及目标区域(interestobject)特性,图层间的权重调配亦有所不同,可依使用的需求加以调整。
2)计算对象的形状(spatial)异质性指标hshape
hshape=ωsmoothness×hsmoothnesscompatness×hcompactness(2)
其中形状的异质性指标是由平滑(smoothness)与紧密(compactness)这两个子异质性指标所构成,ωsmoothness与ωcompactness代表两者间的权重调配,两者的和为1;平滑指标与紧密指标计算如下式所示:
h smoothness = n Merge × l Merge b Merge - ( n Obj 1 × l Obj 1 b Obj 1 + n Obj 2 × l Obj 2 b Obj 2 ) - - - ( 3 )
h compactness = n Merge × l Merge n Merge - ( n Obj 1 × l Obj 1 n Obj 1 + n Obj 2 × l Obj 2 n Obj 2 ) - - - ( 4 )
其中,l为对象的实际边长,b为对象的最短边长,n为对象面积;若平滑指标的权重较高,分割后的对象边界较为平滑,反之,若紧密指标的权重较高,分割后的对象形状较为紧密较接近矩形,根据不同的影像特性以及目标对象(interestobject)特性,两者间的权重调配亦有所不同,可依使用者的需求加以调整。加入形状的因子于影像分割的过程中,能约制对象形状的发展,使分割后的区域形状较平滑完整,较符合人的视觉习惯。
3)计算对象的整体异质性指标h:
h=ωcolor×hcolorshape×hshape(5)
对象整体的异质性指标是由上述的光谱异质性指标与形状异质性指标所构成,ωcolor与ωshape代表光谱与形状两者间的权重调配,两者之和为1,亦可依使用的需求进行调整。对象整体的异质性指标必须满足:
h≤hsc(6)
根据上述的多分辨率分割步骤,分别在尺度l1、在尺度l2、在尺度l3上分别对X1、X2、X1和X2组合的影像进行多尺度分割,分割后的影像分别记为S1、S2和S3
多分辨率分割中几个参数的选择说明如下:
(1)对于高分辨率影像,尺度的选择较大,一般选择50-150。
(2)光谱权重设置为0.8,形状权重设置为0.2。
(3)光滑度权重设置为0.5,紧密度权重设置为0.5。
2a)将S1叠加到X2,根据S1的每一个对象边界确定X2的相对应的对象边界,从而获得X2的分割影像X2S1
2b)将S2叠加到X1,根据S2的每一个对象边界确定X1的相对应的对象边界,从而获得X1的分割影像X1S2
2c)将S3分别叠加到X1和X2,根据S3的每一个对象边界分别确定X1和X2的相对应的对象边界,从而获得X1和X2的分割影像X1S3和X2S3
2d)将S1和S2进行并集操作,获得分割影像S12
2e)将S12分别叠加到X1和X2,根据S12的每一个对象边界分别确定X1和X2的相对应的对象边界,从而获得X1和X2的分割影像X1S12和X2S12
步骤3:根据步骤2a)的结果,对所得到的S1和X2S1进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量X,其中:
Xi=(Mi1,…,Mib,Si1,…,Sib)T(7)
式中,i表示对象的序号,b表示影像的波段数目,Mib表示第b个波段第i个对象的均值,
Sib表示第b个波段第i个对象的方差。
假设Xi是高斯分布,其均值矢量为M,其可设置为X的平均值,X的协方差矩阵为Σ,
Ci=(Xi-M)Σ-1(Xi-M)T(8)
式中Ci是当前对象变化向量卡方(χ2)变换结果值。假设其有2b个自由度,则有:
P ( C i < &chi; 1 - &alpha; 2 ( 2 b ) ) = 1 - &alpha; - - - ( 9 )
如果α取值较小,那么假设检验就有1-α置信度表明Ci的值超过α为潜在的变化对象。
α可以手工调整,建议选择为0.05。根据设定的α得到最终的变化和非变化区域。
步骤4:根据步骤2b)的结果,对所得到的S2和X1S2进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据χ2分布确定影像的变化区域和非变化区域。
步骤5:根据步骤2c)的结果,对所得到的X1S3和X2S3进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据χ2分布确定影像的变化区域和非变化区域。
步骤6:根据步骤2e)的结果,对所得到的X1S12和X2S12进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域。
步骤7:将步骤3-6的检测结果进行基于投票法的融合,获得最终的检测结果。

Claims (3)

1.一种基于融合策略的面向对象的高分辨率遥感影像变化检测方法,其特征在于,包括下述步骤:
步骤1:输入两时相已配准的高分辨率光学遥感影像,分别记为:X1和X2
步骤2:对输入的高分辨率影像进行多分辨率分割:
2a)分别在尺度l1、在尺度l2、在尺度l3上分别对X1、X2、X1和X2组合的影像进行多尺度分割,分割后的影像分别记为S1、S2和S3
2b)将S1叠加到X2,根据S1的每一个对象边界确定X2的相对应的对象边界,从而获得X2的分割影像X2S1
2c)将S2叠加到X1,根据S2的每一个对象边界确定X1的相对应的对象边界,从而获得X1的分割影像X1S2
2d)将S3分别叠加到X1和X2,根据S3的每一个对象边界分别确定X1和X2的相对应的对象边界,从而获得X1和X2的分割影像X1S3和X2S3
2e)将S1和S2进行并集操作,获得分割影像S12
2f)将S12分别叠加到X1和X2,根据S12的每一个对象边界分别确定X1和X2的相对应的对象边界,从而获得X1和X2的分割影像X1S12和X2S12
步骤3:对所得到的S1和X2S1进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域;
步骤4:对所得到的S2和X1S2进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域;
步骤5:对所得到的X1S3和X2S3进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域;
步骤6:对所得到的X1S12和X2S12进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量,并根据卡方分布确定影像的变化区域和非变化区域;
步骤7:将步骤3-6的检测结果进行基于投票法的融合,获得最终的检测结果。
2.如权利要求1所述的基于融合策略的面向对象的高分辨率遥感影像变化检测方法,其特征在于,步骤2中,对输入的高分辨率影像进行多分辨率分割,在尺度hsc下,采用的多分辨率分割包括如下步骤:
1)计算对象的光谱异质性指标hcolor
h c o l o r = &Sigma; c w c &CenterDot; &sigma; c - - - ( 1 )
其中ωc为图层的权重,σc为图层的标准差,c为图层数;根据不同的影像特性以及目标区域特性,图层间的权重调配亦有所不同,依使用的需求加以调整;
2)计算对象的形状异质性指标hshape
hshape=ωsmoothness×hsmoothnesscompatness×hcompactness(2)
其中形状的异质性指标是由平滑与紧密这两个子异质性指标所构成,ωsmoothness与ωcompactness代表两者间的权重调配,两者的和为1;平滑指标与紧密指标计算如下式所示:
h s m o o t h n e s s = n M e r g e &times; l M e r g e b M e r g e - ( n O b j 1 &times; l O b j 1 b O b j 1 + n O b j 2 &times; l O b j 2 b O b j 2 ) - - - ( 3 )
h c o m p a c t n e s s = n M e r g e &times; l M e r g e n M e r g e - ( n O b j 1 &times; l O b j 1 n O b j 1 + n O b j 2 &times; l O b j 2 n O b j 2 ) - - - ( 4 )
其中,l为对象的实际边长,b为对象的最短边长,n为对象面积;若平滑指标的权重较高,分割后的对象边界较为平滑,反之,若紧密指标的权重较高,分割后的对象形状较为紧密较接近矩形,根据不同的影像特性以及目标对象特性,两者间的权重调配亦有所不同,可依使用者的需求加以调整;加入形状的因子于影像分割的过程中,能约制对象形状的发展,使分割后的区域形状较平滑完整,较符合人的视觉习惯;
3)计算对象的整体异质性指标h:
h=ωcolor×hcolorshape×hshape(5)
对象整体的异质性指标是由上述的光谱异质性指标与形状异质性指标所构成,ωcolor与ωshape代表光谱与形状两者间的权重调配,两者之和为1,ωcolor与ωshape依使用的需求进行调整;对象整体的异质性指标必须满足:
h≤hsc(6)
根据上述的多分辨率分割步骤,分别在尺度l1、在尺度l2、在尺度l3上分别对X1、X2、X1和X2组合的影像进行多尺度分割,分割后的影像分别记为S1、S2和S3
3.如权利要求1所述的基于融合策略的面向对象的高分辨率遥感影像变化检测方法,其特征在于,步骤3中,对所得到的S1和X2S1进行相减,并计算每一个对象的均值和方差;将所有对象的均值和方差组合形成一个新的向量X,其中:
Xi=(Mi1,…,Mib,Si1,…,Sib)T(7)
式中,i表示对象的序号,b表示影像的波段数目,Mib表示第b个波段第i个对象的均值,Sib表示第b个波段第i个对象的方差;
假设Xi是高斯分布,其均值矢量为M,其设置为X的平均值,X的协方差矩阵为∑,
Ci=(Xi-M)∑-1(Xi-M)T(8)
式中Ci是当前对象变化向量卡方(χ2)变换结果值;假设其有2b个自由度,则有:
P ( C i < &chi; 1 - &alpha; 2 ( 2 b ) ) = 1 - &alpha; - - - ( 9 )
如果α取值较小,那么假设检验就有1-α置信度表明Ci的值超过α为潜在的变化对象;
其中α取值较小指的是α=0.05;根据设定的α得到最终的变化和非变化区域。
CN201310460801.1A 2013-09-30 2013-09-30 基于融合策略的面向对象的高分辨率遥感影像变化检测方法 Expired - Fee Related CN103489193B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310460801.1A CN103489193B (zh) 2013-09-30 2013-09-30 基于融合策略的面向对象的高分辨率遥感影像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310460801.1A CN103489193B (zh) 2013-09-30 2013-09-30 基于融合策略的面向对象的高分辨率遥感影像变化检测方法

Publications (2)

Publication Number Publication Date
CN103489193A CN103489193A (zh) 2014-01-01
CN103489193B true CN103489193B (zh) 2016-07-06

Family

ID=49829391

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310460801.1A Expired - Fee Related CN103489193B (zh) 2013-09-30 2013-09-30 基于融合策略的面向对象的高分辨率遥感影像变化检测方法

Country Status (1)

Country Link
CN (1) CN103489193B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104036499B (zh) * 2014-05-30 2017-02-22 中国科学院遥感与数字地球研究所 一种多尺度叠加分割方法
CN105160355B (zh) * 2015-08-28 2018-05-15 北京理工大学 一种基于区域相关和视觉单词的遥感图像变化检测方法
CN109934799B (zh) * 2016-09-09 2021-06-25 南京工程学院 多时相差异影像模值计算及变化检测方法
CN106897679A (zh) * 2017-02-13 2017-06-27 长江水利委员会长江科学院 一种基于改进模糊c均值聚类的语义变化检测方法及系统
CN107194942B (zh) * 2017-03-27 2020-11-10 广州地理研究所 一种确定影像分类分割尺度阈值的方法
CN107330875B (zh) * 2017-05-31 2020-04-21 河海大学 基于遥感图像正反向异质性的水体周边环境变化检测方法
CN113361554B (zh) * 2020-03-06 2024-05-24 北京眼神智能科技有限公司 生物特征识别多模态融合方法、装置、存储介质及设备
CN113963222B (zh) * 2021-10-28 2022-09-02 中国电子科技集团公司第五十四研究所 一种基于多策略组合的高分辨率遥感影像变化检测方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102254319A (zh) * 2011-04-19 2011-11-23 中科九度(北京)空间信息技术有限责任公司 一种多层次分割的遥感影像变化检测方法
CN103077525A (zh) * 2013-01-27 2013-05-01 西安电子科技大学 基于Treelet图像融合的遥感图像变化检测方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8620093B2 (en) * 2010-03-15 2013-12-31 The United States Of America As Represented By The Secretary Of The Army Method and system for image registration and change detection

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102254319A (zh) * 2011-04-19 2011-11-23 中科九度(北京)空间信息技术有限责任公司 一种多层次分割的遥感影像变化检测方法
CN103077525A (zh) * 2013-01-27 2013-05-01 西安电子科技大学 基于Treelet图像融合的遥感图像变化检测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A Comparsion of Four Alaorithms for Change Detection in an Urban Environment;Merrill K. Ridd 等;《Remote sensing of Environment》;19981231;95-100 *
基于影像分割的SAR影像变化检测技术研究;王刚;《中国优秀硕士学位论文全文数据库信息科技辑》;20130615;44-60 *

Also Published As

Publication number Publication date
CN103489193A (zh) 2014-01-01

Similar Documents

Publication Publication Date Title
CN103489193B (zh) 基于融合策略的面向对象的高分辨率遥感影像变化检测方法
CN102938066B (zh) 一种基于多元数据重建建筑物外轮廓多边形的方法
CN103839267B (zh) 一种基于形态学建筑物指数的建筑物提取方法
CN105956557B (zh) 一种面向对象的时序遥感影像云覆盖区域自动检测方法
CN106683112B (zh) 一种基于高分辨率图像的道路路域建筑物变化提取方法
CN107844743A (zh) 一种基于多尺度分层残差网络的图像多字幕自动生成方法
CN104915672B (zh) 一种基于高分辨率遥感图像的矩形建筑物提取方法及系统
CN109816641A (zh) 基于多尺度形态学融合的加权局部熵红外小目标检测方法
CN103632363A (zh) 基于多尺度融合的对象级高分辨率遥感影像变化检测方法
CN104240256A (zh) 一种基于层次化稀疏建模的图像显著性检测方法
CN102214291A (zh) 一种快速准确的基于视频序列的人脸检测跟踪方法
CN103761747B (zh) 一种基于加权分布场的目标跟踪方法
CN104103076A (zh) 一种基于高分辨率遥感影像的核电厂规划限制区遥感监查方法
CN104197933A (zh) 一种望远镜视场中高星等滑动恒星的增强和提取方法
Cheng et al. Multi-scale Feature Fusion and Transformer Network for urban green space segmentation from high-resolution remote sensing images
CN105894513A (zh) 顾及影像对象时空变化的遥感影像变化检测方法及系统
Zhang et al. Semantic segmentation of very high-resolution remote sensing image based on multiple band combinations and patchwise scene analysis
CN106373120A (zh) 基于非负矩阵分解和核fcm的多时相遥感影像变化检测方法
Li et al. Detecting building changes using multi-modal Siamese multi-task networks from very high resolution satellite images
CN106384352A (zh) 基于融合策略和fcm的多时相遥感影像变化检测方法
Wang et al. Remote sensing image semantic segmentation method based on small target and edge feature enhancement
CN106021610A (zh) 一种基于显著区域的视频指纹提取方法
CN112487977A (zh) 一种基于国产gf2遥感影像的机场跑道提取方法
Wang Design of public building space in smart city based on big data
Qin et al. Dual-stream infrared object detection network combined with saliency images

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20160706

Termination date: 20170930