CN108804388A - 一种基于eemd的hht的太阳黑子面积周期特征分析的方法 - Google Patents
一种基于eemd的hht的太阳黑子面积周期特征分析的方法 Download PDFInfo
- Publication number
- CN108804388A CN108804388A CN201810478895.8A CN201810478895A CN108804388A CN 108804388 A CN108804388 A CN 108804388A CN 201810478895 A CN201810478895 A CN 201810478895A CN 108804388 A CN108804388 A CN 108804388A
- Authority
- CN
- China
- Prior art keywords
- frequency
- sunspot
- areas
- data
- mode function
- 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.)
- Granted
Links
- 206010064127 Solar lentigo Diseases 0.000 title claims abstract description 29
- 238000000034 method Methods 0.000 title claims abstract description 29
- 230000000737 periodic effect Effects 0.000 title claims abstract description 14
- 238000004458 analytical method Methods 0.000 claims abstract description 22
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 13
- 238000009826 distribution Methods 0.000 claims description 10
- 238000001228 spectrum Methods 0.000 claims description 7
- 239000004744 fabric Substances 0.000 claims description 4
- 238000000926 separation method Methods 0.000 claims description 4
- 230000010355 oscillation Effects 0.000 claims description 2
- 230000003044 adaptive effect Effects 0.000 abstract description 2
- 238000012545 processing Methods 0.000 abstract description 2
- 230000000694 effects Effects 0.000 description 5
- 230000009466 transformation Effects 0.000 description 5
- 230000007774 longterm Effects 0.000 description 3
- 238000012216 screening Methods 0.000 description 3
- 230000014509 gene expression Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
- 230000005428 wave function Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Complex Calculations (AREA)
Abstract
本发明涉及一种基于EEMD的HHT的太阳黑子面积周期特征分析的方法，属天文技术和信号处理领域。本发明采用基于集成经验模式分解(EEMD)算法的希尔伯特黄变换(HHT)的太阳黑子面积周期特征分析的方法，把不同时间尺度上的波动从原始数据中分离出来，不同频率随时间的变化也被清晰的表现出来，与常用的方法相比，它是直接的，自适应的，基于数据本身的方法，其有效地解决了传统的周期分析方法受基函数约束以及时间尺度不准确的问题。
Description
技术领域
本发明涉及一种基于EEMD的HHT的太阳黑子面积周期特征分析的方法，属天文技术和信号处理领域。
背景技术
太阳黑子是太阳活动的一种基本现象，研究太阳黑子的性质，对理解太阳内部结构、空间环境等问题有着重要的意义。人们通常用太阳黑子的相对数或黑子面积来研究它的变化规律。其中，黑子面积数记录每天出现在可见日面的黑子的面积，某种意义上讲，它记录的是一种太阳磁活动的每日磁流量，表征的是太阳发电机产生黑子的功率，以太阳半球面积的百万分之一为单位。相关研究表明，黑子面积的变化体现了太阳日面活动(如光斑、谱斑、日珥、暗条及耀斑等)的变化。和黑子数相比较，黑子面积更具有物理意义，是刻画许多太阳活动的重要物理量。
从黑子面积的观测以来，为了研究黑子面积特性，已经提出很多种分析方法，我们发现这些分析方法主要分析黑子面积的周期特征以及黑子面积的预测，关于黑子面积的周期特征分析方法，例如：(1)傅里叶变换：将时间序列数据从时域变换到频域进行分析研究，太阳活动在不同时间尺度上的变化以其在频域里得到的幅度峰值表现出来。(2)小波分析：基于傅里叶变换，其主要思想是平移和伸缩，可以得到时间序列周期和频率特征的局部精细结构，与傅里叶相比小波分析可以体现变换的详细过程。通过相关文献对这几种方法的介绍，可以看到这些方法对黑子面积的特征分析方法存在着以下问题：(1)傅里叶变换将太阳活动在不同时间尺度上的变化以其在频域里得到的幅度峰值表现出来，并不能得到完全分离的某个时间尺度上的变化分量。而太阳活动在不同时间尺度上变化的周期性对于研究太阳的长期演化具有重要意义。(2)对于长期太阳活动的光谱分析，使用傅里叶变换和其他更先进的基于傅里叶的技术(小波分析)，都有着与基函数先验公式相关联的缺点，例如谐波函数。这样的方法高度限制它们在非平稳信号中的应用。为了解决上述几种方法存在的不足，我们对太阳黑子面积本身进行了分析发现：黑子面积的测量晚于黑子数的观测，1874年开始有连续的黑子面积测量与整理，至今有100余年的观测数据，包括12个完整的太阳活动周及两个半个的太阳活动周。黑子面积数据是长期的非平稳信号。基于这些因素，使得到目前为止各个黑子面积周期特征分析算法都不能很好的分析出黑子周期特性。
发明内容
本发明提供了一种基于EEMD的HHT的太阳黑子面积周期特征分析的方法，以用于解决传统的周期分析方法受基函数约束及时间尺度不准确的问题。
本发明的技术方案是：一种基于EEMD的HHT的太阳黑子面积周期特征分析的方法，所述方法的具体步骤如下：
步骤1、对太阳黑子面积数据进行平滑处理，得到平滑后的数据；
步骤2：将平滑后的数据进行集成经验模式分解，分解后得到从高频到低频的一系列不同频率的本质模态函数集合；其中本质模态函数集合包括所有的本质模态函数分量及残差，每个本质模态函数分量代表着太阳黑子面积信号中的一个固有的振动模态；
步骤3：根据时间间隔确定采样频率，根据数据长度和采样频率得到频率分辨率，由奈奎斯特定理得出最大、最小频率值，对本质模态函数集合中的每一个IMF分量进行希尔伯特黄变换，得到其对应的瞬时频率谱；
步骤4：统计每一个瞬时频率谱中每个频率出现的次数，将统计的频率出现的次数的最大值变为1，并以该值为基准对所有频率的出现次数作归一化处理，得到对应的频率概率分布直方图；
步骤5：用高斯函数对频率概率分布直方图进行单高斯拟合，拟合出的高斯曲线的峰值的横坐标作为该本质模态函数的平均频率，取一个标准差的置信区间，得到该本质模态函数的频率区间；
步骤6：步骤5中得到的平均频率的倒数为平均周期，平均频率加一个标准差的和的倒数为该本质模态函数的最小周期，平均频率减去一个标准差的差的倒数为该本质模态函数的最大周期。
所述步骤2中，集成经验模式分解参数：添加噪声等级Nstd＝0.2，添加噪声次数NR＝300，最大本质模态函数个数numIMF＝7，最大迭代筛选次数maxSift＝10，得到IMF1-IMF6和残差。
所述步骤3中，希尔伯特黄变换参数：采样频率fs＝1/150，频率分辨率f0＝fs/N，由奈奎斯特定理，最大频率fmax＝fs/2，最小频率fmin＝f0；其中，N为平滑后数据的长度。
本发明的有益效果是：本发明采用基于集成经验模式分解(EEMD)算法的希尔伯特黄变换(HHT)的太阳黑子面积周期特征分析的方法，把不同时间尺度上的波动(即，本质模态函数，IMF)从原始数据中分离出来，不同频率随时间的变化也被清晰的表现出来，与常用的方法相比，它是直接的，自适应的，基于数据本身的方法，其有效地解决了传统的周期分析方法受基函数约束以及时间尺度不准确的问题。
附图说明
图1是本发明基于集成经验模式分解算法的希尔伯特黄变换的太阳黑子面积周期特征分析的方法的总体流程图；
图2是本发明中采用1874年5月1日至2016年10月31日(日平均)的南半球太阳黑子面积数据；
图3是本发明中对图2经平滑处理后的图像；
图4是本发明中对图3进行集成经验模式分解得到的本质模态函数集合(IMFs)与残差(R)；
图5、6是本发明中图4中IMF1分量的瞬时频率和高斯拟合后对应的频率概率分布梯度图；
图7、8是本发明中图4中IMF2分量的瞬时频率和高斯拟合后对应的频率概率分布梯度图；
图9、10是本发明中图4中IMF3分量的瞬时频率和高斯拟合后对应的频率概率分布梯度图；
图11、12是本发明中图4中IMF4分量的瞬时频率和高斯拟合后对应的频率概率分布梯度图；
图13、14是本发明中图4中IMF5分量的瞬时频率和高斯拟合后对应的频率概率分布梯度图；
图15、16是本发明中图4中IMF6分量的瞬时频率和高斯拟合后对应的频率概率分布梯度图。
具体实施方式
实施例1：如图1-16所示，一种基于EEMD的HHT的太阳黑子面积周期特征分析的方法，所述方法的具体步骤如下：
步骤1：数据平滑：首先对太阳黑子面积数据进行平滑处理(采用1874年5月至2016年10月(日平均)的南半球太阳黑子面积数据；图2为南半球太阳黑子面积的原始数据，图3为平滑后的数据图像)。
SouthSmooth＝smooth(South,w,'moving')
其中，SouthSmooth是平滑后的数据，South是太阳黑子原始数据，w为平滑尺度，moving表示用移动平均滤波器的方法来平滑数据，n为第n个数据，将连续的采样数据看成一个长度固定为N的队列，在新的一次测量后，上述队列的首数据去掉，其余N-1个数据依次前移，并将新的采样数据插入，作为新队列的尾；然后对这个队列进行算术运算。本发明中w＝150。
步骤2：使用集成经验模式分解方法分解数据：
(1)将平滑后的数据添加噪声，改善信号的极值点分布，通过借助高斯白噪声的频率均布特性，在频域为信号构建频率均布的尺度。添加噪声公式如下：
x(t)＝s(t)+n(t)
其中，x(t)为添加噪声后的数据，s(t)为真实信号，n(t)为添加的噪声。
(2)将添加完噪声的数据分解为IMFs。这种分解过程又被称作“筛选”。首先，找出添加完噪声的数据的所有局部极值点(极大值，极小值)，分别对极大值点、极小值点进行三次样条插值函数拟合形成数据的上下包络线，求出上下包络线的均值，再用原始数据减去上下包络线的均值得到初始IMF，具体公式如下：
x(t)-m1＝h1
其中，s(x)为三次样条插值函数，hi＝xi+1-xi，i＝0,1,...,n-1，Mi＝s”(xi)为待定参数，M0，M1…Mn满足线性方程。插值后得到的上包络线为emax，下包络线为emin，上下包络的平均值为m1，s(x)为三次样条插值函数，x(t)同上，初始的IMF为h1。
(3)使用不同的白噪声，重复(1)、(2)步骤。过程如下：
h1-m11＝h11
h1(k-1)-m1k＝h1k
fmax＝fs/2
对h1数据进行接下来的一轮筛选，重复筛选过程k次，如果h1k满足IMF分量的两个条件，那么得到第一个本质模态函数分量c1。IMF分量的两个判断条件为：1.在整个数据域内，局部极值点与过零点数目必须相等或至多相差一个；2.在任意点，由局部极值构成的上下包络线的平均值为零。
(4)对应的IMFs集合的均值作为分解的最终结果。
分解得到的结果满足：
这里x(t)同上，n为分解出来的IMF分量的个数，cj表示第j个IMF分量，rn表示残差。得到的IMFs集合与残差，如图4所示。
所述步骤2中，集成经验模式分解参数：添加噪声等级Nstd＝0.2，添加噪声次数NR＝300，最大本质模态函数个数numIMF＝7，最大迭代筛选次数maxSift＝10，得到IMF1-IMF6和残差。
步骤3：对每一个本质模态函数进行希尔伯特变换：首先根据时间间隔确定采样频率，根据数据长度和采样频率得到频率分辨率，由奈奎斯特定理得出最大、最小频率值，由希尔伯特黄变换得到瞬时频率谱。具体公式如下：
[nt,tscale,fscale]＝nnspa(IMFs(i,:),t0,t1,fres,tres,fmin,fmax,tw0,tw1,'hilbert','none')
其中，nt为光谱的二维矩阵表示，tscale为表示时间轴值的向量，fscale表示频率轴值的向量，IMFs(i,:)表示对第i个IMF分量进行希尔伯特黄变换，t0为开始时间，t1为结束时间，fres表示频率分辨率，tres为时间分辨率，fmin,fmax为满足奈奎斯特定理的最大最小频率，tw0,tw1表示时间缩放尺度，这里不设置缩放，tw0,tw1与对应的t0t1相等。hilbert表示使用希尔伯特变换，none表示不标准化。希尔伯特变换的结果为Y(t)，其中P为柯西主值，x(t)为给定的时间序列(同上)。本发明中t0＝0，t1＝N/fs，fres＝fmax/fmin，tres＝N，fmin＝fs/N，fmax＝fs/2，fs＝1/w，本发明中w＝150，平滑后的数据长度N＝347。
步骤4：统计每一个瞬时频率谱中每个频率出现的次数，将统计的频率出现的次数的最大值变为1，并以该值为基准对所有频率的出现次数作归一化处理(如，假设有频率1出现100次；频率2出现200次；则将200变为1，100变为0.5)，得到对应的频率概率分布直方图；
步骤5：使用高斯拟合频率概率分布直方图：用高斯函数对频率分布概率直方图进行单高斯拟合，拟合出的高斯曲线的峰值的横坐标作为该本质模态函数的平均频率，取一个标准差的置信区间[b-c b+c]。得到该本质模态函数的频率区间。具体公式如下：y(x)＝aexp([x-b]2/2c2)；其中待估参数a、b、c分别为高斯曲线的峰值、峰值对应的横坐标及标准差。也就是说，a为平均频率出现的次数，b为平均频率，c为标准差。
步骤3、4、5的结果，如图5至图16所示。举例说明IMF1分量的情况，其它等同；其中图6中的梯度线为图5瞬时频率的概率统计，虚线曲线为高斯拟合曲线，高斯拟合曲线峰值代表的点分别对应频率轴、概率轴拟合的平均频率和平均频率出现的概率。垂直于频率轴(频率轴为过峰值且垂直于横坐标的轴)的两条点线为一个标准差的置信区间。
步骤6：求取频率的倒数得到周期区间：平均频率的倒数为平均周期，平均频率加一个标准差的和的倒数为该本质模态函数的最小周期，平均频率减去一个标准差的差的倒数为该本质模态函数的最大周期。有平均周期为1/b，周期置信区间为[1/(b+c)1/(b-c)]。最终的周期统计结果，如表1所示(表1是本发明中图4中所有IMFs分量对应的平均频率、平均周期(天、年不同单位)、及周期置信区间)。
上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。
Claims (3)
1.一种基于EEMD的HHT的太阳黑子面积周期特征分析的方法，其特征在于：所述方法的具体步骤如下：
步骤1、对太阳黑子面积数据进行平滑处理，得到平滑后的数据；
步骤2：将平滑后的数据进行集成经验模式分解，分解后得到从高频到低频的一系列不同频率的本质模态函数集合；其中本质模态函数集合包括所有的本质模态函数分量及残差，每个本质模态函数分量代表着太阳黑子面积信号中的一个固有的振动模态；
步骤3：根据时间间隔确定采样频率，根据数据长度和采样频率得到频率分辨率，由奈奎斯特定理得出最大、最小频率值，对本质模态函数集合中的每一个IMF分量进行希尔伯特黄变换，得到其对应的瞬时频率谱；
步骤4：统计每一个瞬时频率谱中每个频率出现的次数，将统计的频率出现的次数的最大值变为1，并以该值为基准对所有频率的出现次数作归一化处理，得到对应的频率概率分布直方图；
步骤5：用高斯函数对频率概率分布直方图进行单高斯拟合，拟合出的高斯曲线的峰值的横坐标作为该本质模态函数的平均频率，取一个标准差的置信区间，得到该本质模态函数的频率区间；
步骤6：步骤5中得到的平均频率的倒数为平均周期，平均频率加一个标准差的和的倒数为该本质模态函数的最小周期，平均频率减去一个标准差的差的倒数为该本质模态函数的最大周期。
2.根据权利要求1所述的基于EEMD的HHT的太阳黑子面积周期特征分析的方法，其特征在于：所述步骤2中，集成经验模式分解参数：添加噪声等级Nstd＝0.2，添加噪声次数NR＝300，最大本质模态函数个数numIMF＝7，最大迭代筛选次数maxSift＝10，得到IMF1-IMF6和残差。
3.根据权利要求1所述的基于EEMD的HHT的太阳黑子面积周期特征分析的方法，其特征在于：所述步骤3中，希尔伯特黄变换参数：采样频率fs＝1/150，频率分辨率f0＝fs/N，由奈奎斯特定理，最大频率fmax＝fs/2，最小频率fmin＝f0；其中，N为平滑后数据的长度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810478895.8A CN108804388B (zh) | 2018-05-18 | 2018-05-18 | 一种基于eemd的hht的太阳黑子面积周期特征分析的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810478895.8A CN108804388B (zh) | 2018-05-18 | 2018-05-18 | 一种基于eemd的hht的太阳黑子面积周期特征分析的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108804388A true CN108804388A (zh) | 2018-11-13 |
| CN108804388B CN108804388B (zh) | 2021-12-17 |
Family
ID=64092650
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810478895.8A Active CN108804388B (zh) | 2018-05-18 | 2018-05-18 | 一种基于eemd的hht的太阳黑子面积周期特征分析的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108804388B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110851627A (zh) * | 2019-09-24 | 2020-02-28 | 昆明理工大学 | 一种用于描述全日面图像中太阳黑子群的方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101527698A (zh) * | 2009-04-03 | 2009-09-09 | 北京理工大学 | 基于希尔伯特黄变换和自适应陷波的非平稳干扰抑制方法 |
| CN105510711A (zh) * | 2015-12-24 | 2016-04-20 | 合肥工业大学 | 一种改进的经验模态分解的谐波分析法 |
| US20170311870A1 (en) * | 2014-11-14 | 2017-11-02 | Neurochip Corporation, C/O Zbx Corporation | Method and apparatus for processing electroencephalogram (eeg) signals |
-
2018
- 2018-05-18 CN CN201810478895.8A patent/CN108804388B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101527698A (zh) * | 2009-04-03 | 2009-09-09 | 北京理工大学 | 基于希尔伯特黄变换和自适应陷波的非平稳干扰抑制方法 |
| US20170311870A1 (en) * | 2014-11-14 | 2017-11-02 | Neurochip Corporation, C/O Zbx Corporation | Method and apparatus for processing electroencephalogram (eeg) signals |
| CN105510711A (zh) * | 2015-12-24 | 2016-04-20 | 合肥工业大学 | 一种改进的经验模态分解的谐波分析法 |
Non-Patent Citations (3)
| Title |
|---|
| NORDENE.HUANG等著: "《希尔伯特黄变换及其应用 第2版》", 31 October 2017, 国防工业出版社 * |
| 刘爽 等: "太阳黑子相对数平滑月均值多尺度周期分析", 《中国海洋大学学报》 * |
| 独舞的深海精灵: "希尔伯特黄变换_百度百科2017-05-24历史版本", 《HTTPS://BAIKE.BAIDU.COM/HISTORY/希尔伯特黄变换/8286826/118809838》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110851627A (zh) * | 2019-09-24 | 2020-02-28 | 昆明理工大学 | 一种用于描述全日面图像中太阳黑子群的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108804388B (zh) | 2021-12-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Hong | Decomposition and forecast for financial time series with high-frequency based on empirical mode decomposition | |
| Kolláth et al. | Multiple and changing cycles of active stars-I. Methods of analysis and application to the solar cycles | |
| CN104751000B (zh) | 一种机电复合传动状态监测信号小波降噪方法 | |
| CN113325277A (zh) | 一种局部放电处理方法 | |
| Benetazzo et al. | A/D converter performance analysis by a frequency-domain approach | |
| Wang et al. | A reducing iteration orthogonal matching pursuit algorithm for compressive sensing | |
| CN103236825A (zh) | 一种用于高精度数据采集系统的数据校正方法 | |
| CN116865269A (zh) | 一种风电机组高谐波补偿方法及系统 | |
| CN110632386B (zh) | 一种太阳射电干扰滤除方法、可读存储介质及电子设备 | |
| CN108804388A (zh) | 一种基于eemd的hht的太阳黑子面积周期特征分析的方法 | |
| CN113189634B (zh) | 一种类高斯成形方法 | |
| González-Rodríguez et al. | Computational method for extracting and modeling periodicities in time series | |
| Das et al. | Review of adaptive decomposition-based data preprocessing for renewable generation rich power system applications | |
| CN116046968A (zh) | 一种液相色谱工作站数据处理方法、系统及可存储介质 | |
| CN111008356B (zh) | 一种基于WTSVD算法扣除背景的γ能谱集分析方法 | |
| Reidy | An introduction to random processes for the spectral analysis of speech data | |
| Ortega et al. | Hilbert–Huang transform analysis of storm waves | |
| CN115128661A (zh) | 一种适用于解析计算γ能谱的方法 | |
| CN110117798B (zh) | 一种铝电解的氧化铝浓度估计方法及装置 | |
| CN112325868A (zh) | 一种基于多尺度变换的偏振光罗盘去噪方法 | |
| Shi et al. | Recurrent neural network wind power prediction based on variational modal decomposition improvement | |
| Kuang et al. | Reconstructing signal from quantized signal based on singular spectral analysis | |
| Li et al. | Bi-dimensional variational mode decomposition for surface texture analysis | |
| Ansoult | Circular sampling for Fourier analysis of digital terrain data | |
| CN109030452A (zh) | 一种基于五点三次平滑算法的拉曼光谱数据降噪方法 |
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 |