CN105116442B - 岩性油气藏弱反射地震信号的重构方法 - Google Patents

岩性油气藏弱反射地震信号的重构方法 Download PDF

Info

Publication number
CN105116442B
CN105116442B CN201510439378.6A CN201510439378A CN105116442B CN 105116442 B CN105116442 B CN 105116442B CN 201510439378 A CN201510439378 A CN 201510439378A CN 105116442 B CN105116442 B CN 105116442B
Authority
CN
China
Prior art keywords
signal
vibration component
seismic
natural mode
instantaneous
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
CN201510439378.6A
Other languages
English (en)
Other versions
CN105116442A (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.)
Yangtze University
Original Assignee
Yangtze University
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 Yangtze University filed Critical Yangtze University
Priority to CN201510439378.6A priority Critical patent/CN105116442B/zh
Publication of CN105116442A publication Critical patent/CN105116442A/zh
Application granted granted Critical
Publication of CN105116442B publication Critical patent/CN105116442B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开一种岩性油气藏弱反射地震信号的重构方法,通过对原始地震信号分别进行小波去噪、波形匹配端点延拓处理,并运用经验模态分解对原始地震信号进行分解得到固有模态分量,采取相关性自适应选择有效固有模态分量做希尔伯特变换,对得到的瞬时振幅进行自适应谱白化处理,并与希尔伯特变换得到的瞬时相位重构信号,得到含噪少能量高的重构地震信号,其相比于原始地震信号不仅提高了信噪比和地震剖面的分辨率,更方便检测与识别。

Description

岩性油气藏弱反射地震信号的重构方法
技术领域
本发明涉及一种地震信号的处理方法,具体涉及一种岩性油气藏弱反射地震信号的重构方法。
背景技术
随着油气能源需求的增加,油气开发的重心逐步由构造油气藏转向岩性油气藏,这对岩性油气藏的地震数据处理技术提出了更高的要求。与常见的构造油气藏相比,岩性油气藏大多都储量分散、层薄、埋藏深、孔渗条件复杂,采集到地震数据具有信号反射弱、含噪大的特点。强反射地震信号容易被处理识别,弱反射地震信号很难处理与识别。
目前,针对弱反射地震信号的检测与识别,已有的处理方法效果都不是很理想。传统的弱反射地震信号处理方法中,傅立叶变换、短时傅立叶变换都受Heisenberg测不准原理制约,不能同时达到时间和频率的高精度要求,故这些变换对非线性非平稳的弱反射地震信号不太适用。后期发展起来的小波变换虽然能够去除弱反射地震信号中的大部分噪声,但是小波变换本身存在缺陷,使得去噪后的弱反射地震信号同相轴边缘变得模糊不清、分辨率下降,在表达二维信号时会损失边缘细节信息,造成剖面成像模糊,分辨率降低。并且小波变换不具备自适应的特点,制约了它在处理弱反射地震信号领域的发展。
1998年,黄锷及其合作者提出了希尔伯特-黄变换(Hilbert-Huang Transform)这一新的针对非线性、非平稳信号进行时频分析的有效方法。希尔伯特-黄变换根据信号本身的特性进行经验模态分解,对分解得到的固有模态分量做希尔伯特变换,能够得到信号有效的瞬时振幅、瞬时相位、瞬时频率。希尔伯特-黄变换不受Heisenberg测不准原理制约,它能在时间和频率同时达到很高的精度,这使得它适于分析弱反射地震信号。近年来,有国内学者提出将希尔伯特-黄变换运用到处理地震信号的领域来,并取得了一些进展。油气开采的重心由构造油气藏转向岩性油气藏,但构造油气藏的处理方法直接在岩性油气藏运用却没有取得较好的效果,阻碍了岩性油气藏的勘探开采。由此可见,现有技术中仍存在诸多缺失,急待加以改良。
发明内容
有鉴于此,有必要提供一种信噪比高、且地震剖面的分辨率良好的地震信号处理方法。
一种岩性油气藏弱反射地震信号的重构方法,所述弱反射地震信号重构方法包括以下步骤:
S1、对原始地震信号进行小波变换的预处理去噪;
S2、将去噪后的信号进行波形匹配端点延拓处理;
S3、对端点延拓处理后的地震信号进行经验模态分解,得到能够反映各种信号特征的固有模态分量;
S4、对分解得到的固有模态分量与原始地震信号进行相关性计算,将得到的相关系数与阈值相关系数进行对比,自适应选择得到有效固有模态分量;
S5、对有效固有模态分量进行希尔伯特变换得到瞬时振幅、瞬时相位、瞬时频率,将瞬时振幅进行白化处理,并对白化后的瞬时振幅和希尔伯特变换得到的瞬时相位重构信号。
本发明所述的岩性油气藏弱反射地震信号的重构方法,通过对原始地震信号分别进行小波去噪、波形匹配端点延拓处理,并运用经验模态分解对原始地震信号进行分解得到固有模态分量,采取相关性自适应选择有效固有模态分量做希尔伯特变换,对得到的瞬时振幅进行自适应谱白化处理,并与希尔伯特变换得到的瞬时相位重构信号,得到含噪少能量高的重构地震信号,其相比于原始地震信号不仅提高了信噪比和地震剖面的分辨率,更方便检测与识别。
附图说明
图1为本发明所述岩性油气藏弱反射地震信号的重构方法的流程框图;
图2为一组原始地震剖面与重构地震剖面对比图;
图3为另一组原始地震剖面与重构地震剖面对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明,应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
本发明所述一种岩性油气藏弱反射地震信号的重构方法,所述弱反射地震信号重构方法包括以下步骤:
S1、对原始地震信号进行小波变换的预处理去噪;
S2、将去噪后的信号进行波形匹配端点延拓处理;
S3、对端点延拓处理后的地震信号进行经验模态分解,得到能够反映各种信号特征的固有模态分量;
S4、对分解得到的固有模态分量与原始地震信号进行相关性计算,将得到的相关系数与阈值相关系数进行对比,自适应选择得到有效固有模态分量;
S5、对有效固有模态分量进行希尔伯特变换得到瞬时振幅、瞬时相位、瞬时频率,将瞬时振幅进行白化处理,并对白化后的瞬时振幅和希尔伯特变换得到的瞬时相位重构信号。
下面采用本发明所述岩性油气藏弱反射地震信号的重构方法对某工区实际地震数据进行处理,以此对本发明作具体的说明与验证。其中,本发明是在Matlab上进行运行操作。具体的,如表1所示:
表1实验测试环境
所述工区的地震数据进为叠后地震数据,共7条剖面,每条剖面地震信号道数为100。本地震数据为岩性油气藏地震勘探弱反射信号,信号受随机噪声干扰较大,存在部分同相轴难识别、地层剖面信息模糊不清、分辨率低的特点。
步骤S1:读取原始地震数据文件中的实际地震信号到计算机内存,对原始地震信号进行小波变换的预处理去噪。
具体的,首先调用地震数据文件操作函数,先读取原始地震数据的文件头,然后读取文件中数据部分,得到原始地震道数、采样点数、采样间隔、炮点坐标等全部信息,其中,每个地震道由240个字节的道头和地震数据组成;设读取的原始地震信号记为h1(t)。
在Matlab中,将原始地震信号的数据采用矩阵操纵,对原始地震信号进行小波去噪的预处理,去除原始地震信号中的随机噪声。小波去噪后得到的信号记为h2(t)。
同时还可以在该步骤中计算出原始地震信号的信噪比,及显示剖面,可用来与重构后的地震信号进行对比。
由于地震信号采集技术的提高,面波等规律噪声已在信号采集过程中得到较好的处理,故不做考虑。但由于弱反射地震信号有用信号能量较弱,随机噪声可能掩盖有用信号的现象会对后期相关性的选取产生影响,故采取小波去噪法对数据做一个预处理。
步骤S2:对于上述经验模态分解过程中会遇到的端点失真的问题,本方法对步骤S1中进行小波去噪后的地震信号利用波形匹配的端点延拓方法进行延拓信号处理,用于提高后续步骤中的分解结果的准确度。以左端点为例,右端点延拓的方法与左端点相同,具体的,所述步骤S2包括以下分步骤:
首先,计算波形匹配度:α(d0,dN,PN),其中,d0是信号中左端第一个最小值点到起始端的波形,其长度记为l,dN是信号中极大值点到左边起始位置的波形,PN是信号中的极大值点。
然后,令α(d0,dN0,PN0)=min{α(d0,dN,PN),N=0,1,2...}。若α(d0,dN0,PN0)<β°l满足,其中β为一常数,则用从dN0起左端的第一个极大值和第一个极小值的信号作为原始信号的左端延拓;若不满足上述条件,则取原始信号左端起始处的第一个和第二个极大值点的均值,作为延拓信号的极大值;取左端起始处的第一个和第二个极小值点的均值,作为延拓信号的极小值,进行延拓。右端延拓方法与左端相同。
步骤S3:由于信号之间的固有特征存在差异,故利用经验模态分解的方法,自适应的将信号分解成若干个能反映各种信号特征的固有模态分量。因此,对上步进行延拓后信号进行经验模态分解,具体步骤如下:
初始化:对信号和信号的分解次数进行编号,设信号与分解次数的编号分别为k和i,令s1(t)=h3(t),i=1,k=1;
②构造待筛分信号:k=k+1,ψk(t)=Sk(t);
③找出ψk(t)中所有的极大值和极小值,用三次样条插值函数求上下包络线,计算上下包络线均值μk(t);
④令ψk+1(t)=ψk(t)-μk(t);
⑤计算阈值限定值判断sd≤0.25是否成立,如果成立,则有固有模态分量,如不成立,则返回步骤③,令ψk(t)=ψk+1(t);
⑥保存第i阶固有模态分量:ci(t)=ψk(t),sk+1=sk(t)-ci(t),i=i+1;
⑦判断Sk+1(t)是否单调或为常量,如何是,则余量记为分解剩余项rk,分解过程结束,否则返回步骤②,令Sk+1(t)=Sk(t)。
经过经验模态分解后的信号可以表示为:
其中h4(t)为被分解后的信号;ci为分解得到的固有模态分量;a为分解产生固有模态的个数;rk为分解剩余项。
步骤S4:由于分解得到的固有模态分量存在伪分量与噪声集中分量,信号理论上分解出a项基本模式分量,但是由于分解误差会产生a1项基本模式分量c″i与a2项伪分量c″i,即:
经验模态分解后的基本模式分量c′i与分解前信号的相关性R如下:
在公式(4)中,由于经验模态分解过程是局部正交分解,所以:
伪分量ck与分解前信号的相关性R为:
由公式(4)(5)(6)可知,各基本模式分量与分解前信号的相关性约等于各分量的自相关,而伪分量与分解前信号的相关性约等于0,所以,可以利用此相关性判别真伪分量达到自适应选择去除伪分量。根据随机噪声的无方向性可知随机噪声集中的分量与原始信号的相关性较低,趋近于0。根据上面的推论我们就可以做一个规定:选取相关系数大于0.01的固有模态分量作为有效分量,得到k个有效固有模态分量Xk(t)。
因为经验模态分解过程中会分出伪分量或噪声集中的固有模态分量,故这些分量需要除去,在此,本方法建立固有模态分量与原始地震信号相关系数大小对比的基础上,来判别各个固有模态分量是否是所需的有用信号集中的分量。
步骤S5:固有模态分量是平稳信号或简单的非线性信号,属于窄带信号。任何一个窄带信号X(t)都能得到它的希尔伯特变换结果Y(t),其公式为:
对步骤S4中的有效固有模态分量Xk(t),进行希尔伯特变换得Yk(t)。
Xk(t)与Yk(t)构成解析信号Zk(t)
其中ak(t)为瞬时振幅,θk(t)为瞬时相位,ωk(t)为瞬时频率。
设计白化滤波器,对有效固有模态分量的瞬时振幅进行谱白化处理,实现瞬时振幅的自适应增强。
首先通过傅里叶变换得到有效固有模态分量的瞬时振幅谱Ak(ω),并求出振幅谱的包络λk(ω),则白化滤波器为:
其中mk为包络λk(ω)的最大值,即mk=max(λk(ω))。ε为白噪因子,用于调节地震数据的信噪比和分辨率。ε越大,分辨率降低,信噪比增强;ε越小,分辨率提高,信噪比降低。ε的数值是根据实际的地震信号来计算选择的。在岩性油气藏弱反射地震信号中,常选择的白噪因子ε数值为0.3。
然后将每个有效固有模态分量上的瞬时振幅ak(t)与白化滤波器fk(ω)相乘,其中的ω由瞬时频率ωk(t)替代,从而实现瞬时振幅的自适应增强。白化后的瞬时振幅为:
对瞬时振幅进行自适应增强后,地震信号的相位信息均保持不变。可以将瞬时相位与增强后的瞬时振幅进行重构,得到重构地震信号。
进一步的,将重构地震信号的信噪比与原始地震信号的信噪比进行比较,其比较结果如表2所示:
表2原始地震剖面与重构地震剖面信噪比
由表2可以看出,重构地震信号的信噪比相较原始地震信号的信噪比提高了,这就说明原始地震中的噪声得到有效压制。
同时,如图2和图3所示,重构地震信号的剖面与原始地震信号的剖面相对比,重构后地震信号的剖面中同向轴不连续的地方已经变得清楚,整个剖面的分辨率得到了提高,同时还提高了储层识别能力,因此,重构的地震信号相较原始地震信号含噪少能量高。
本发明所述的岩性油气藏弱反射地震信号的重构方法,通过对原始地震信号分别进行小波去噪、波形匹配端点延拓处理,并运用经验模态分解对原始地震信号进行分解得到固有模态分量,采取相关性自适应选择有效固有模态分量做希尔伯特变换,对得到的瞬时振幅进行自适应谱白化处理,并与希尔伯特变换得到的瞬时相位重构信号,得到含噪少能量高的重构地震信号,其相比于原始地震信号不仅提高了信噪比和地震剖面的分辨率,更方便检测与识别。
以上所述仅为本发明的较佳实施例,并不用以限制本发明,凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (5)

1.一种岩性油气藏弱反射地震信号的重构方法,所述弱反射地震信号重构方法包括以下步骤:
S1、对原始地震信号进行小波变换的预处理去噪;
S2、将去噪后的信号进行波形匹配端点延拓处理;
所述步骤S2包括以下子步骤:
S21、计算波形匹配度:α(d0,dN,PN);
其中,d0是信号中左端第一个最小值点到起始端的波形,其长度记为l,dN是信号中极大值点到左边起始位置的波形,PN是信号中的极大值点,令α(d0,dN0,PN0)=min{α(d0,dN,PN),N=0,1,2...};
S22、若α(d0,dN0,PN0)<β°l满足,其中β为一常数,则用从dN0起左端的第一个极大值和第一个极小值的信号作为原始地震信号的左端延拓;
S23、若不满足上述条件,则取原始信号左端起始处的第一个和第二个极大值点的均值,作为延拓信号的极大值;取左端起始处的第一个和第二个极小值点的均值,作为延拓信号的极小值,进行延拓;
S3、对端点延拓处理后的地震信号进行经验模态分解,得到能够反映各种信号特征的固有模态分量;
S4、对分解得到的固有模态分量与原始地震信号进行相关性计算,将得到的相关系数与阈值相关系数进行对比,自适应选择得到有效固有模态分量;
S5、对有效固有模态分量进行希尔伯特变换得到瞬时振幅、瞬时相位、瞬时频率,将瞬时振幅进行白化处理,并对白化后的瞬时振幅和希尔伯特变换得到的瞬时相位重构信号。
2.根据权利要求1所述的一种岩性油气藏弱反射地震信号的重构方法,其特征在于:所述步骤S1包括以下子步骤:
S11、读取地震记录数据的文件头,并读取地震记录数据文件中的数据部分,得到地震信号的各项信息;
S12、将数据进行矩阵操纵,对原始地震信号进行小波去噪。
3.根据权利要求1所述的一种岩性油气藏弱反射地震信号的重构方法,其特征在于:所述步骤S3包括以下子步骤:
S31、将延拓后信号初始化,找出延拓后信号的所有的极大值和极小值点;
S32、通过三次样条插值函数求上下包络线,计算上下包络的均值,通过包络线的均值得到筛选门限的阈值限定值;
S33、判断阈值限定值是否不处于设定的阈值区间范围内,若不在阈值区间范围内,则有固有模态分量,否则,返回步骤S31,对下一个信号进行循环判断;
S34、判断余量是否为单调函数或是常量,如果是,则余量记为分解剩余项,如果不是,则返回步骤S31,对下一个信号进行循环判断。
4.根据权利要求1所述的一种岩性油气藏弱反射地震信号的重构方法,其特征在于:所述步骤S4包括以下子步骤:
S41、设置有效固有模态分量的阈值相关系数,将经验模态分解后得到的固有模态分量与分解前信号进行相关性比对;
S42、如果比对结果大于阈值相关系数,则该固有模态分量为有效固有模态分量;
S43、如果比对结果小于阈值相关系数,则该固有模态分量为伪分量。
5.根据权利要求1所述的一种岩性油气藏弱反射地震信号的重构方法,其特征在于:所述步骤S5包括以下子步骤:
S51、对筛选出的有效固有模态分量进行希尔伯特变换,变换得到的瞬时振幅、瞬时相位、瞬时频率;
S52、通过傅里叶变换得到有效固有模态分量的瞬时振幅谱,并求出振幅谱的包络,求得白化滤波器;
S53、将每个有效固有模态分量的瞬时振幅与白化滤波器相乘,实现瞬时振幅的自适应增强;
S54、通过步骤S51中的瞬时相位与S53中自适应增强后的瞬时振幅进行重构信号。
CN201510439378.6A 2015-07-24 2015-07-24 岩性油气藏弱反射地震信号的重构方法 Expired - Fee Related CN105116442B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510439378.6A CN105116442B (zh) 2015-07-24 2015-07-24 岩性油气藏弱反射地震信号的重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510439378.6A CN105116442B (zh) 2015-07-24 2015-07-24 岩性油气藏弱反射地震信号的重构方法

Publications (2)

Publication Number Publication Date
CN105116442A CN105116442A (zh) 2015-12-02
CN105116442B true CN105116442B (zh) 2019-01-01

Family

ID=54664474

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510439378.6A Expired - Fee Related CN105116442B (zh) 2015-07-24 2015-07-24 岩性油气藏弱反射地震信号的重构方法

Country Status (1)

Country Link
CN (1) CN105116442B (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106371138B (zh) * 2016-08-17 2018-10-16 中国石油天然气集团公司 地震数据重建方法和装置
CN106405444B (zh) * 2016-08-29 2019-01-29 东南大学 一种基于改进的经验模态分解的地磁测量信号去噪方法
CN108072894A (zh) * 2016-11-18 2018-05-25 中国石油化工股份有限公司 一种地震信号处理方法
CN107515424A (zh) * 2017-07-26 2017-12-26 山东科技大学 一种基于vmd与小波包的微震信号降噪滤波方法
CN108051857B (zh) * 2017-12-05 2019-08-06 中国石油大学(华东) 一种基于经验模态分解的局部相似性小波阈值去噪方法
CN109164490A (zh) * 2018-09-10 2019-01-08 中国石油化工股份有限公司 一种提高沙漠地区弱反射成像精度的方法
CN109613609B (zh) * 2019-01-16 2019-09-24 国家深海基地管理中心 一种基于变长度处理的组合信号分解方法
CN112200229A (zh) * 2020-03-02 2021-01-08 周红 基于传感器通信的信号处理方法及装置
CN111551992B (zh) * 2020-05-12 2021-02-26 中国科学院地质与地球物理研究所 岩石储层构造表征方法、装置、计算机可读存储介质及电子设备
CN112200069B (zh) * 2020-09-30 2022-11-04 山东大学 时频域谱减法和经验模态分解联合的隧道滤波方法及系统
CN112378414A (zh) * 2020-11-20 2021-02-19 深圳信息职业技术学院 一种基于pm2.5健康出行的路径规划装置及方法
CN113608018A (zh) * 2021-06-30 2021-11-05 中冶南方都市环保工程技术股份有限公司 提高谐波检测精度的自适应vmd检测方法、装置及存储介质

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102998706B (zh) * 2012-11-23 2015-05-20 中国石油大学(北京) 一种衰减地震数据随机噪声的方法及系统
CN103308940B (zh) * 2013-05-27 2016-08-10 中山大学 地震剖面的经验模态分解方法
CN104182617B (zh) * 2014-07-30 2017-04-12 电子科技大学 基于内禀波形匹配的端点效应抑制方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
一种自适应的EMD端点延拓方法;邵晨曦 等;《电子学报》;20071031;第35卷(第10期);第1944-1948页
利用改进希尔伯特—黄变换进行地震资料时频分析;曹思远 等;《石油地球物理勘探》;20130430;第48卷(第2期);第246-254页

Also Published As

Publication number Publication date
CN105116442A (zh) 2015-12-02

Similar Documents

Publication Publication Date Title
CN105116442B (zh) 岩性油气藏弱反射地震信号的重构方法
Gómez et al. A simple method inspired by empirical mode decomposition for denoising seismic data
Battista et al. Application of the empirical mode decomposition and Hilbert-Huang transform to seismic reflection data
Wu et al. Noise attenuation for 2-D seismic data by radial-trace time-frequency peak filtering
CN107102356B (zh) 基于ceemd的地震信号高分辨率处理方法
CN107356964B (zh) S变换域基于变分原理的q值估计与补偿方法
CN106199532B (zh) 基于混合傅立叶-小波分析的探地雷达信号降噪方法
CN109633753B (zh) 地表噪声压制方法及装置
CN108020863A (zh) 一种基于地震奇偶函数的碳酸盐岩薄储层孔隙度预测方法
CN107179550B (zh) 一种数据驱动的地震信号零相位反褶积方法
CN106934183B (zh) 频散曲线确定方法和装置,及纵横波速度确定方法和装置
Alsalmi et al. Mask filtering to the Wigner-Ville distribution
CN112183407B (zh) 一种基于时频域谱减法的隧道地震波数据去噪方法及系统
Han et al. Firefly algorithm with disturbance-factor-based particle filter for seismic random noise attenuation
CN106291733B (zh) 一种便携式石油检测器
Li et al. Removing abnormal environmental noise in nodal land seismic data using deep learning
Hao et al. Denoising method based on spectral subtraction in time-frequency domain
CN110727025A (zh) 一种隐蔽断层的识别方法
CN108646296A (zh) 基于自适应谱峭度滤波器的沙漠地震信号噪声消减方法
Guili et al. Decomposition of rock micro-fracture signals based on a singular value empirical mode decomposition algorithm
CN112713907B (zh) 一种基于字典学习的海洋csem降噪方法及系统
CN112764108B (zh) 一种基于改进经验小波变换的新型地震资料噪声压制算法
Xu et al. Signal filtering based on wavelet transform and its application in ground penetrating radar
CN110568491B (zh) 一种品质因子q的估算方法
CN107730494A (zh) 一种基于变分模态分解的锚杆检测方法

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20190101

Termination date: 20210724

CF01 Termination of patent right due to non-payment of annual fee