CN109829872B - 一种用于内陆水体遥感的多时相多源遥感影像融合方法 - Google Patents

一种用于内陆水体遥感的多时相多源遥感影像融合方法 Download PDF

Info

Publication number
CN109829872B
CN109829872B CN201910060135.XA CN201910060135A CN109829872B CN 109829872 B CN109829872 B CN 109829872B CN 201910060135 A CN201910060135 A CN 201910060135A CN 109829872 B CN109829872 B CN 109829872B
Authority
CN
China
Prior art keywords
image
time
images
abundance
spectrum
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
CN201910060135.XA
Other languages
English (en)
Other versions
CN109829872A (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.)
Henan Agricultural University
Original Assignee
Henan Agricultural 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 Henan Agricultural University filed Critical Henan Agricultural University
Priority to CN201910060135.XA priority Critical patent/CN109829872B/zh
Publication of CN109829872A publication Critical patent/CN109829872A/zh
Application granted granted Critical
Publication of CN109829872B publication Critical patent/CN109829872B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Processing (AREA)

Abstract

本发明涉及用于内陆水体遥感的多时相多源遥感影像融合方法，有效解决同时融合多源影像的空间、光谱、时间分辨率，生成适用于内陆光学复杂水体遥感监测的数据源，用于对水体的监测和环保的问题，利用高空间、光谱分辨率影像，以及高时间分辨率影像，获得影像数据，并进行预处理，生成高光谱、高时间分辨率影像，混合光谱分解，得到高时间频率的丰度图；抽取低时频影像成像时刻的丰度数据，结合高空间、光谱分辨率影像，得到高光谱分辨率的端元光谱，基于高时频丰度图和端元光谱，生成初步融合结果，将GOCI影像和初步融合结果重新排列，利用初步融合结果及高空间、光谱分辨率影像生成最终融合结果，本发明方法易操作，应用面广，精度高，稳定可靠。

Description

一种用于内陆水体遥感的多时相多源遥感影像融合方法
技术领域
本发明涉及遥感技术，特别是一种用于内陆水体遥感的多时相多源遥感影像融合方法。
背景技术
遥感技术地质勘探、测量中常用的技术，因其能够快速、准确的对水体状况进行检测，具有传统定点观测无法比拟的优势，因此已经成为水环境监测的重要手段。对内陆光学复杂水体而言，首先，其水体边界复杂，需要影像具有较高的空间分辨率；其次，水体光学性质复杂，欲对其进行较精准的监测，需要影像具有较高的光谱分辨率；最后，由于水体状况变化迅速，因此需要影像具有较高的时间分辨率。受限于遥感传感器硬件条件，目前没有传感器能同时具备以上特性。理想数据源的缺失限制了对内陆光学复杂水体的遥感监测。
多源影像融合技术为该问题提供了一种可行的方案。目前的图像融合算法可分为三种：第一种，基于色彩变换、以及小波变换等数学变换的空间-光谱融合算法。通过对这些图像融合方法进行分类、比较和评估发现，这些算法只能对准同步多源图像进行融合。也就是说，空间-光谱融合算法丢失了一些时间数据(信息)。第二种，利用STARFM算法将LandsatTM和MODIS反射率数据进行空间-时间融合，并用于区域土地覆盖制图，变化监测，作物生物量的估算，这种算法用于将高时间、低空间分辨率图像与低时间、高空间分辨率图像(如具有LandsatTM图像的MODIS图像)融合。从光谱上看，只有具有相似波长的波段才能被融合。也就是说，空间-时间融合算法丢失了一些光谱数据(信息)。第三种，空间-光谱-时间融合算法，基于空间细节数据的时间变化可忽略的假设，可以融合多源影像的空间、光谱以及时间分辨率。但首先，这类算法中空间细节数据可以忽略的假设对于水体来说显然不成立，其次，这种算法能无法覆盖所有情况：例如地球静止卫星图像融合(GOCI和Himawari图像)，其中空间和光谱数据都由低时间分辨率图像提供。因此，影像融合上的改进和创新势在必行。
发明内容
针对上述情况，为克服现有技术之缺陷，本发明之目的就是提供一种用于内陆水体遥感的多时相多源遥感影像融合方法，可有效解决同时融合多源影像的空间、光谱、时间分辨率，以生成适用于内陆光学复杂水体遥感监测的数据源，有效用于对水体的监测和环保的问题
本发明解决的技术方案是，一种用于内陆水体遥感的多时相多源遥感影像融合方法，包括以下步骤：
(1)利用高空间、光谱分辨率影像，以及高时间分辨率影像，获得影像数据，并进行预处理，生成高光谱、高时间分辨率影像，首先利用SMACC算法将高时间分辨率影像进行混合光谱分解，得到高时间频率的丰度图；之后抽取低时频影像成像时刻的丰度数据，结合高空间、光谱分辨率影像，利用下式计算得到高光谱分辨率的端元光谱：
Figure BDA0001953868740000021
其中，Rλ(l,s)是像元(l,s)处，波长λ处的反射率值，C(l,s；m)是第m个端元的丰度图，
Figure BDA0001953868740000022
是第m个端元的反射率值，S和L分别代表图像的行列数，J代表分度求解目标函数；
2)、初步融合：基于高时频丰度图和高光谱分辨率的端元光谱，生成初步融合结果，初步融合结果具有2km空间分辨率，8个光谱波段以及10min的时间分辨率；
分两步对GOCI、Himawari影像进行融合，首先，利用SMACC算法将Himawari影像进行混合光谱分分解，得到Himawari端元光谱和丰度图像；之后抽取8:30-15:30之间，以1h为间隔的丰度图，利用该丰度图和GOCI影像，通过最小二乘法最小化中的J，计算得到GOCI端元光谱：
Figure BDA0001953868740000023
其中，Rλ(l,s)是像元(l,s)处，波长λ处的反射率值，C(l,s；m)是第m个端元的丰度图，
Figure BDA0001953868740000024
是第m个端元的反射率值，S和L分别代表图像的行列数；
基于Himawari丰度图和GOCI端元，通过式(4)计算得到初步融合结果：
Figure BDA0001953868740000025
其中，
Figure BDA0001953868740000026
代表估算得到的反射率值，
Figure BDA0001953868740000027
是计算得到的丰度图，
Figure BDA0001953868740000028
是计算得到的端元光谱；
3)、将GOCI影像和初步融合结果进行重新排列，得到时域影像，利用初步融合结果以及高空间、光谱分辨率影像生成最终融合结果：首先将空域影像转化到时域，将原图像逐波段叠加的形式打散，将多时相的数据进行叠加生成时域影像，时域影像每个波段代表一个时段的影像数据，不同波段的数据被分为不同的影像；在时域中，利用SMACC算法那将初步融合结果影像进行分解，得到高时频端元光谱以及丰度图；再基于高空间、光谱分辨率影像重采样得到的高光谱分辨率影像和上一步得到的端元数据(信息)，利用最小二乘法计算得到高时频端元光谱；随后，利用上一步得到的端元数据(信息)以及高空间、光谱分辨率影像，利用下式计算得到高空间分辨率的丰度图：
Figure BDA0001953868740000031
其中，F代表端元求解目标函数，将高空间分辨率的丰度图与高时频的端元数据(信息)逐像元相乘，利用下式计算得到最终融合结果：
Figure BDA0001953868740000032
其中，
Figure BDA0001953868740000033
代表估算得到的反射率值，
Figure BDA0001953868740000034
是计算得到的丰度图，
Figure BDA0001953868740000035
是计算得到的端元光谱。
本发明方法简单，易操作，应用面广，精度高，稳定可靠，对图像的随机噪声以及椒盐噪声都有较好的控制效果，可广泛应用于内陆水环境遥感监测，在当前传感器硬件条件受限的情况下为水环境遥感提供更加理想的数据源，具有很强的实际应用价值，非常有利于内陆水体的监测以及环境保护，有显著的经济和社会效益
附图说明
图1为本发明方法的框示图。
具体实施方式
以下结合附图和具体情况对本发明的具体实施方式作详细说明。
本发明在具体实施中，一种用于内陆水体遥感的多时相多源遥感影像融合方法，包括以下步骤：
1)、利用高空间、光谱分辨率影像，以及高时间分辨率影像，获取影像数据，并进行预处理：
对获取的影像数据利用图像自带的经纬度网格进行几何校正，将图像的投影坐标都转化为UTM投影，WGS-84坐标系，将影像获取的信号转化为大气层外反射率，GOCI和Himawari影像分别利用式(1)和式(2)进行转换：
Figure BDA0001953868740000036
其中，RTOA是大气顶层反射率，LTOA是GOCI影像获取得到的大气层外下行辐亮度，F0为大气层外的太阳辐照度，λ是波长，θ是太阳高度角，在GOCI影像中，太阳高度角以图层的方式存储；
Figure BDA0001953868740000041
其中，albedo(λ)是Himawari影像获取到的大气层外二向反射率数据；
2)、初步融合，即GOCI和Himawari影像融合：
分两步对GOCI、Himawari影像进行融合，首先，利用SMACC算法将Himawari影像进行混合光谱分分解，得到Himawari端元光谱和丰度图像；之后抽取8:30-15:30之间，以1h为间隔的丰度图，利用该丰度图和GOCI影像，通过最小二乘法最小化式(3)中的J，计算得到GOCI端元光谱：
Figure BDA0001953868740000042
其中，Rλ(l,s)是像元(l,s)处，波长λ处的反射率值，C(l,s；m)是第m个端元的丰度图，
Figure BDA0001953868740000043
是第m个端元的反射率值，S和L分别代表图像的行列数；
基于Himawari丰度图和GOCI端元，通过式(4)计算得到初步融合结果：
Figure BDA0001953868740000044
其中，
Figure BDA0001953868740000045
代表估算得到的反射率值，
Figure BDA0001953868740000046
是计算得到的丰度图，
Figure BDA0001953868740000047
是计算得到的端元光谱；
3)、将GOCI影像和初步融合结果进行重排列，得到时域影像；利用SMACC算法将初步融合结果影像进行分解，得到初步融合端元光谱以及丰度图；进而，基于GOCI影像重采样得到2kmGOCI图像，结合上一步得到的端元数据(信息)，利用式(3)得到高时频端元光谱；随后，利用上一步得到的高时频端元光谱，以及GOCI影像，利用最小二乘法最小化式(5)中的F得到500m分辨率丰度图：
Figure BDA0001953868740000048
最后利用式(4)将500m空间分辨率的丰度图与高时频的端元数据(信息)结合，得到最终融合结果。
本发明经实地试验和应用，效果非常好，以内陆太湖水域为例，首先，通过网站http://kosc.kiost.ac/eng/获取GOCI影像数据(如图1)，通过网站http://www.eorc.jaxa.jp/ptree获取Himawari影像数据(如图1)。数据地点是太湖，数据时间是2016年2月22日；并同现有其他几种方法进行了试验对比，以下是本发明的图像融合方法与几种常用方法的对比结果。虽然目前针对GOCI和Himawari这种类型融合的算法还未见报道，但在本发明第二步中，初步融合结果与GOCI的融合，可以看作是高光谱数据与多光谱数据的融合，针对这种融合，已有诸多成熟算法。因此选取多种方法将初步融合结果与GOCI影像进行融合之后，与本发明的方法进行对比。采用GS、NND、PC、UBF等四种方法与初步融合结果(s1)相结合，形成s1GS、s1NND、s1PC、s1UBF四种方法，采用相关系数(r)、平均相对误差(MAPE)以及尺度无关误差(ERGAS)来对融合结果进行评价。参考图像为GOCI影像。
Figure BDA0001953868740000051
其中，Sh是高空间分辨率影像像元大小，Sl是地空间分辨率影像像元大小，Z代表波段个数，i代表当前参与计算的波段，RMSE代表均方根误差，M代表参考影像的均值。ERGAS值越小，代表两景影像越接近。
Figure BDA0001953868740000052
其中，x、y分别是梁静影像的像元值，n是像元数，X、Y分别为均值，r越大代表两景影像相关性越高。
Figure BDA0001953868740000053
MAPE越小代表两景影像的差距越小。
Figure BDA0001953868740000054
从三个精度指标看来，s1GS、s1NND、s1PC以及s1UBF方法相较于融合结果都有一定提升，但本发明方法的ERGAS值和MAPE值均明显低于其他几种方法，改进效果比较明显。说明本发明方法求得较好的效果。由上表可以看出，本发明方法明显优于现有方法，而且相关指标都得到了大幅度改善，表明本发明方法效果非常好，而且经过反复多次试验，都取得了相同或相近似的结果，表明方法稳定、可靠，具有很强的实际应用价值。
以上表明，本发明的多时相多源遥感影像融合方法，以混合光谱分解为基础，在谱域和时域两个方向对多源影像进行分解及重构，最终完成图像融合，融合结果同时具有较高的空间分辨率、光谱分辨率以及时间分辨率。相较于目前常用的融合方法，具有高精度，对误差不敏感等特点。综上，本发明特别适合于内陆光学复杂水体多源遥感影像的融合。对内陆水环境具有重要实际应用价值，在传感器硬件条件难以满足需求的现状下生成更加适用的数据源，可有效用于水体质量的监测和环境保护，经济和社会效益显著。

Claims (2)

1.一种用于内陆水体遥感的多时相多源遥感影像融合方法，其特征在于，包括以下步骤：
(1)利用高空间、光谱分辨率影像，以及高时间分辨率影像，获得影像数据，并进行预处理，生成高光谱、高时间分辨率影像，首先利用SMACC算法将高时间分辨率影像进行混合光谱分解，得到高时间频率的丰度图；之后抽取低时频影像成像时刻的丰度数据，结合高空间、光谱分辨率影像，利用下式计算得到高光谱分辨率的端元光谱：
Figure FDA0001953868730000011
其中，Rλ(l,s)是像元(l,s)处，波长λ处的反射率值，C(l,s；m)是第m个端元的丰度图，
Figure FDA0001953868730000012
是第m个端元的反射率值，S和L分别代表图像的行列数，J代表分度求解目标函数；
2)、初步融合：基于高时频丰度图和高光谱分辨率的端元光谱，生成初步融合结果，初步融合结果具有2km空间分辨率，8个光谱波段以及10min的时间分辨率；
分两步对GOCI、Himawari影像进行融合，首先，利用SMACC算法将Himawari影像进行混合光谱分分解，得到Himawari端元光谱和丰度图像；之后抽取8:30-15:30之间，以1h为间隔的丰度图，利用该丰度图和GOCI影像，通过最小二乘法最小化中的J，计算得到GOCI端元光谱：
Figure FDA0001953868730000013
其中，Rλ(l,s)是像元(l,s)处，波长λ处的反射率值，C(l,s；m)是第m个端元的丰度图，
Figure FDA0001953868730000014
是第m个端元的反射率值，S和L分别代表图像的行列数；
基于Himawari丰度图和GOCI端元，通过式(4)计算得到初步融合结果：
Figure FDA0001953868730000015
其中，
Figure FDA0001953868730000016
代表估算得到的反射率值，
Figure FDA0001953868730000017
是计算得到的丰度图，
Figure FDA0001953868730000018
是计算得到的端元光谱；
3)、将GOCI影像和初步融合结果进行重新排列，得到时域影像，利用初步融合结果以及高空间、光谱分辨率影像生成最终融合结果：首先将空域影像转化到时域，将原图像逐波段叠加的形式打散，将多时相的数据进行叠加生成时域影像，时域影像每个波段代表一个时段的影像数据，不同波段的数据被分为不同的影像；在时域中，利用SMACC算法那将初步融合结果影像进行分解，得到高时频端元光谱以及丰度图；再基于高空间、光谱分辨率影像重采样得到的高光谱分辨率影像和上一步得到的端元数据，利用最小二乘法计算得到高时频端元光谱；随后，利用上一步得到的端元数据以及高空间、光谱分辨率影像，利用下式计算得到高空间分辨率的丰度图：
Figure FDA0001953868730000021
其中，F代表端元求解目标函数，将高空间分辨率的丰度图与高时频的端元数据逐像元相乘，利用下式计算得到最终融合结果：
Figure FDA0001953868730000022
其中，
Figure FDA0001953868730000023
代表估算得到的反射率值，
Figure FDA0001953868730000024
是计算得到的丰度图，
Figure FDA0001953868730000025
是计算得到的端元光谱。
2.根据权利要求1所述的用于内陆水体遥感的多时相多源遥感影像融合方法，其特征在于，包括以下步骤：
1)、利用高空间、光谱分辨率影像，以及高时间分辨率影像，获取影像数据，并进行预处理：
对获取的影像数据利用图像自带的经纬度网格进行几何校正，将图像的投影坐标都转化为UTM投影，WGS-84坐标系，将影像获取的信号转化为大气层外反射率，GOCI和Himawari影像分别利用式(1)和式(2)进行转换：
Figure FDA0001953868730000026
其中，RTOA是大气顶层反射率，LTOA是GOCI影像获取得到的大气层外下行辐亮度，F0为大气层外的太阳辐照度，λ是波长，θ是太阳高度角，在GOCI影像中，太阳高度角以图层的方式存储；
Figure FDA0001953868730000027
其中，albedo(λ)是Himawari影像获取到的大气层外二向反射率数据；
2)、初步融合，即GOCI和Himawari影像融合：
分两步对GOCI、Himawari影像进行融合，首先，利用SMACC算法将Himawari影像进行混合光谱分分解，得到Himawari端元光谱和丰度图像；之后抽取8:30-15:30之间，以1h为间隔的丰度图，利用该丰度图和GOCI影像，通过最小二乘法最小化式(3)中的J，计算得到GOCI端元光谱：
Figure FDA0001953868730000031
其中，Rλ(l,s)是像元(l,s)处，波长λ处的反射率值，C(l,s；m)是第m个端元的丰度图，
Figure FDA0001953868730000032
是第m个端元的反射率值，S和L分别代表图像的行列数；
基于Himawari丰度图和GOCI端元，通过式(4)计算得到初步融合结果：
Figure FDA0001953868730000033
其中，
Figure FDA0001953868730000034
代表估算得到的反射率值，
Figure FDA0001953868730000035
是计算得到的丰度图，
Figure FDA0001953868730000036
是计算得到的端元光谱；
3)、将GOCI影像和初步融合结果进行重排列，得到时域影像；利用SMACC算法将初步融合结果影像进行分解，得到初步融合端元光谱以及丰度图；进而，基于GOCI影像重采样得到2kmGOCI图像，结合上一步得到的端元数据，利用式(3)得到高时频端元光谱；随后，利用上一步得到的高时频端元光谱，以及GOCI影像，利用最小二乘法最小化式(5)中的F得到500m分辨率丰度图：
Figure FDA0001953868730000037
最后利用式(4)将500m空间分辨率的丰度图与高时频的端元数据结合，得到最终融合结果。
CN201910060135.XA 2019-01-22 2019-01-22 一种用于内陆水体遥感的多时相多源遥感影像融合方法 Active CN109829872B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910060135.XA CN109829872B (zh) 2019-01-22 2019-01-22 一种用于内陆水体遥感的多时相多源遥感影像融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910060135.XA CN109829872B (zh) 2019-01-22 2019-01-22 一种用于内陆水体遥感的多时相多源遥感影像融合方法

Publications (2)

Publication Number Publication Date
CN109829872A CN109829872A (zh) 2019-05-31
CN109829872B true CN109829872B (zh) 2022-11-22

Family

ID=66861880

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910060135.XA Active CN109829872B (zh) 2019-01-22 2019-01-22 一种用于内陆水体遥感的多时相多源遥感影像融合方法

Country Status (1)

Country Link
CN (1) CN109829872B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111721714B (zh) * 2020-06-01 2021-07-09 东北林业大学 一种基于多源光学遥感数据的土壤含水量估算方法
CN113160100A (zh) * 2021-04-02 2021-07-23 深圳市规划国土房产信息中心(深圳市空间地理信息中心) 一种基于光谱信息影像的融合方法、融合装置及介质
CN113222836B (zh) * 2021-04-25 2024-03-01 自然资源部国土卫星遥感应用中心 一种高光谱与多光谱遥感信息融合方法及系统
CN116310881B (zh) * 2023-04-24 2023-08-11 中国科学院空天信息创新研究院 土壤有机质含量估算方法、装置、电子设备及存储介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1642087A1 (en) * 2003-07-09 2006-04-05 Fabrizio Ferrucci Method and apparatus for automatically detecting and mapping, particularly for burnt areas without vegetation
CN104809691A (zh) * 2015-05-05 2015-07-29 李云梅 一种基于滑动窗口混合像元分解的图像融合方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1642087A1 (en) * 2003-07-09 2006-04-05 Fabrizio Ferrucci Method and apparatus for automatically detecting and mapping, particularly for burnt areas without vegetation
CN104809691A (zh) * 2015-05-05 2015-07-29 李云梅 一种基于滑动窗口混合像元分解的图像融合方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于生物光学模型的水体多源遥感图像融合算法研究;郭宇龙等;《光学学报》;20150410(第04期);全文 *

Also Published As

Publication number Publication date
CN109829872A (zh) 2019-05-31

Similar Documents

Publication Publication Date Title
CN109829872B (zh) 一种用于内陆水体遥感的多时相多源遥感影像融合方法
Lanaras et al. Super-resolution of Sentinel-2 images: Learning a globally applicable deep neural network
JP7004844B2 (ja) 経験的直交関数分解法による静止海色衛星データの再構築方法
Duggin et al. Assumptions implicit in remote sensing data acquisition and analysis
CN109993237B (zh) 基于高分卫星光学遥感数据的水体快速提取方法及系统
Rakwatin et al. Restoration of Aqua MODIS band 6 using histogram matching and local least squares fitting
Fasbender et al. Bayesian data fusion for adaptable image pansharpening
CN104154998B (zh) 一种基于压缩感知的计算多光谱成像图谱的重构方法
CN204788661U (zh) 基于压缩感知的计算多光谱成像系统
CN108427964B (zh) 一种遥感图像与地球化学的融合方法及系统
Kim et al. Correction of stray-light-driven interslot radiometric discrepancy (ISRD) present in radiometric products of geostationary ocean color imager (GOCI)
Duysak et al. Investigation of the performance of different wavelet-based fusions of SAR and optical images using Sentinel-1 and Sentinel-2 datasets
CN113744249A (zh) 一种海洋生态环境损害调查方法
Guo et al. A novel multitemporal image-fusion algorithm: Method and application to GOCI and himawari images for inland water remote sensing
Jing et al. Sub-pixel accuracy evaluation of FY-3D MERSI-2 geolocation based on OLI reference imagery
Aggarwal et al. Multi-spectral demosaicing technique for single-sensor imaging
Fryskowska et al. Some aspects of satellite imagery integration from Eros B and Landsat 8
Neigh et al. Quantifying Libya-4 surface reflectance heterogeneity with WorldView-1, 2 and EO-1 Hyperion
Fayad et al. Vision Transformers, a new approach for high-resolution and large-scale mapping of canopy heights
Lee et al. Bathymetry of shallow coastal regions derived from space-borne hyperspectral sensor
Kolluru et al. Haze correction of hyperspectral imagery over inland waters
Alvera-Azcárate et al. Creation of high resolution suspended particulate matter data in the North Sea from Sentinel-2 and Sentinel-3 data
Oehmcke et al. Estimating Forest Canopy Height With Multi-Spectral and Multi-Temporal Imagery Using Deep Learning
Liu Bathymetry and bottom albedo retrieval using Hyperion: a case study of Thitu Island and reef
CN116721243B (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