CN111623749B - 一种基于D-InSAR的铁路采空区边界提取新方法 - Google Patents

一种基于D-InSAR的铁路采空区边界提取新方法 Download PDF

Info

Publication number
CN111623749B
CN111623749B CN202010423214.5A CN202010423214A CN111623749B CN 111623749 B CN111623749 B CN 111623749B CN 202010423214 A CN202010423214 A CN 202010423214A CN 111623749 B CN111623749 B CN 111623749B
Authority
CN
China
Prior art keywords
goaf
railway
image
boundary
insar
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
CN202010423214.5A
Other languages
English (en)
Other versions
CN111623749A (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.)
China Railway Design Corp
Original Assignee
China Railway Design Corp
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 China Railway Design Corp filed Critical China Railway Design Corp
Priority to CN202010423214.5A priority Critical patent/CN111623749B/zh
Publication of CN111623749A publication Critical patent/CN111623749A/zh
Application granted granted Critical
Publication of CN111623749B publication Critical patent/CN111623749B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C5/00Measuring height; Measuring distances transverse to line of sight; Levelling between separated points; Surveyors' levels
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021SAR image post-processing techniques
    • G01S13/9023SAR image post-processing techniques combined with interferometric techniques
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/03Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于D‑InSAR的铁路采空区边界提取新方法,包括以下步骤:获取研究区域内时间序列SAR影像和精密轨道文件;将SAR影像数据预处理建立连接图;对连接图中的干涉对进行干涉处理生成差分干涉图,将差分干涉图进行地理编码;利用U‑net神经网络对地理编码后差分干涉图进行采空区识别生成采空区‑非采空区二值图像;对识别后采空区二值图像栅格转矢量;对多期矢量化后的采空区结果进行空间融合;统计分析采空区面积大小、发展方向和采空区边界距铁路的最近距离。本发明对威胁铁路安全的地下采空区进行了自动识别,并对采空区的分布特点和发展趋势进行了系统分析,为优化铁路选线和排查铁路采空区隐患点提供了技术支持。

Description

一种基于D-InSAR的铁路采空区边界提取新方法
技术领域
本发明涉及铁路采空区边界提取新方法,具体的是一种基于D-InSAR(Differential Interferometric Synthetic Aperture Radar差分干涉测量)的铁路采空区边界提取新方法。
背景技术
采空区是一个突出的生态环境问题,它会导致周围产生不同程度的地面沉降,直接对地表的房屋、道路和基础设施的安全造成威胁。在优化铁路走向方案时,需充分考虑安全性和后期维护成本,主动避开这些具有隐患的采空区。因而,快速、准确识别采空区对优化铁路选线方案具有重要意义。目前,进行采空区的监测识别主要采用水准测量、GPS等传统技术手段,虽然可以获得实时、高精度的点位监测结果,但其成果只能是以离散点的形式,且获取周期长,无法快速获取直观的空间立体沉降分布信息,给铁路设计方案的宏观决策带来困难。与之相比,近年兴起的D-INSAR技术,可利用SAR影像全天时、覆盖范围广、穿透能力强的特点,快速获得大面积的沉陷盆地范围和形变场,在矿区沉陷识别监测中得到了很好的应用。然而采空区往往位于地质环境复杂、沉降量级大、变化速度快的矿区,D-INSAR技术无法直接获取大梯度的地形形变信息,其应用受到了限制。
发明内容
因此,本发明的目的在于提供一种基于D-InSAR的铁路采空区边界提取新方法,解决上述背景技术中的技术问题。
为了实现上述目的,本发明的一种基于D-InSAR(Differential InterferometricSynthetic Aperture Radar差分干涉测量)的铁路采空区边界提取新方法,包括以下步骤:
(1)获取研究区域内时间序列SAR影像和精密轨道文件;,
(2)将SAR影像数据预处理(聚焦、配准、裁剪)建立干涉对连接图;
(3)对连接图中的干涉对进行干涉处理生成差分干涉图;
(4)将差分干涉图进行地理编码;
(5)利用U-net神经网络对地理编码后差分干涉图进行采空区识别生成采空区-非采空区二值图像;
(6)对识别后采空区-非采空区二值图像由栅格转矢量,生成矢量化的采空区结果;
(7)对多期矢量化后的采空区结果进行空间融合;
(8)统计分析采空区面积大小、发展方向和采空区边界距铁路的最近距离。
获取研究区域内时间序列SAR影像和精密轨道文件,具体步骤为:每月时间段内所获取的SAR影像不少于2幅同极化SLC影像。同时,获取对应SAR影像的卫星精密轨道参数文件。
将SAR影像数据预处理建立连接图,具体步骤为:将所有SAR影像聚焦、配准到超级主影像,对影像进行去斜,用精密轨道文件优化影像参数文件中的轨道参数,根据矢量文件裁剪覆盖研究区的范围,再对所有影像建立干涉对连接图。
影像空间垂直基线为临界基线的0-30%,时间基线为0-95d,每幅影像干涉对数应至少保证2对。
将未自动参与组网影像按时空最邻近原则人工添加参与组网。
所述的干涉处理具体步骤为:对建立连接图的干涉对进行差分干涉、去平、轨道精炼、重去平、滤波处理;所述的地理编码差分干涉图是将差分干涉图从雷达坐标系转换到WGS84地理坐标系。
利用U-net神经网络对地理编码后图像进行采空区识别生成采空区-非采空区二值图像,具体步骤为:利用U-net神经网络结构,对编码后的差分干涉图进行语义分割,对采空区图像进行像素级别的分类,从而进行采空区识别生成采空区-非采空区二值图像。
对识别后采空区-非采空区二值图像栅格转矢量,具体步骤为:将U-net提取的二分类栅格采空区-非采空区二值图像,根据二分类图像的梯度变化,将栅格采空区图像进行矢量化处理。
所述的对多期矢量化后的采空区结果进行空间融合,具体步骤为:针对多期矢量化后的采空区结果,用空间拓扑分析进行矢量融合,将每期的小范围矢量取并集融合成整体采空区矢量。
统计分析采空区面积大小、发展方向和采空区边界距铁路的最近距离,具体步骤为:对融合后的整体采空区矢量进行统计分析,获取采空区面积大小,发展方向以及采空区边缘距铁路的最近距离。
采用上述技术方案,本发明的有益技术效果为:
本发明基于D-INSAR技术特点,可以利用Sentinel-1A/B卫星影像求取差分干涉图,通过U-net神经网络从差分干涉图上自动提取铁路采空区,经矢量融合获取采空区分布特点、面积大小、发展方向以及与铁路的最近距离信息,对威胁铁路安全的地下采空区进行了自动识别,并对采空区的分布特点和发展趋势进行了系统分析,为优化铁路选线和排查铁路采空区隐患点提供了技术支持。
附图说明
图1为矢量融合示意图。
图2为一种基于D-InSAR的铁路采空区边界提取新方法的流程图。
图3为SAR影像干涉组网连接图。
具体实施方式
以下通过附图和具体实施方式对本发明作进一步的详细说明。
如图所示,本发明的一种结合附图对本发明的技术方案作进一步说明。如图2所示,本发明所涉及的一种基于D-InSAR的铁路采空区边界提取新方法包括以下技术步骤:
S1、获取研究区域内时间序列SAR影像和精密轨道文件。
Sentinel-1A/B卫星搭载C波段合成孔径雷达,是目前现势较好的免费SAR数据,其影像分辨率高,双星重访周期为6天,极其适合采空区探测。考虑到采空区地质条件复杂、沉降量级大、变化速度快等特点,因而所获取的雷达影像在时间上连续性要好,每月时间段内所获取的SAR影像不少于2幅同极化SLC影像,影像选择同极化方式的单视复数SLC影像。同时,下载对应SAR影像日期对应的卫星精密轨道参数文件。
S2、数据预处理建立干涉对连接图。
干涉处理前将SAR影像数据预处理建立连接图,将所有SAR影像聚焦、配准到超级主影像,根据矢量文件裁剪覆盖研究区的影像,对影像进行去斜,用精密轨道文件优化影像参数文件中的轨道参数。再对所有SAR影像进行干涉对组网,建立干涉对连接图(见图3);考虑到时间周期长相干性差,生成的干涉图条纹信息不明显、图像噪声明显的问题,组网应尽量符合短空间基线和短时间基线的原则,保证所有影像空间垂直基线为临界基线的0-30%,时间基线为0-95d,此外每幅影像干涉对数应至少保证2对。在上述原则下,若有少数影像未参与组网,通过人工手动添加参与组网,建立最终的连接图文件。
S3、对干涉对连接图进行干涉处理生成差分干涉图,将差分干涉图进行地理编码。
所述的干涉处理具体步骤为:对建立连接图的影像进行差分干涉、去平、轨道精炼、重去平、滤波处理;经过S2处理建立所有影像连接图,根据连接图文件进行干涉处理。干涉处理需要去除平地效应,引入外部SRTM 30m分辨率的DEM,建立查询列表将WGS84坐标的DEM地理编码到SAR影像所在的雷达坐标系下,DEM采样因子控制在2。轨道精炼通过在地形平坦、没有相位跃变和形变条纹的区域均匀选取GCP控制点,采用三次多项式估算轨道相位偏移量,同时基于控制点对干涉相位进行再去平。为了消除干涉图斑点噪声,采用自适应滤波对干涉图进行滤波处理。在处理过程中,根据数据类型需要设置流程处理的相关参数,其中Azimuth Looks设置为1,Range Looks设置为4,Cartographic Grid Size设置为20m,重采样方法为四次立体卷积法,滤波因子为0.65,滤波窗口为64。得到的差分干涉图是基于雷达坐标系的,为了便于后续沉陷区矢量属性信息的提取与分析,需要将差分干涉图从雷达坐标系转换到WGS84地理坐标系。借助干涉处理中建立的查询列表和DEM参数文件,将差分干涉图反向编码,从雷达坐标系转换到WGS84地理坐标系。
S4、利用U-net神经网络对编码后图像进行采空区识别生成采空区-非采空区二值图像。
利用U-net神经网络对地理编码后图像进行采空区识别生成采空区图像,具体步骤为:利用U-net神经网络结构,对编码后的差分干涉图进行语义分割,对采空区图像进行像素级别的分类,从而进行采空区识别生成采空区-非采空区二值图像。
U-net网络基于Encoder-Decoder结构,结构简明且稳定,所需训练样本少。
在进行模型训练前,对S3地理编码后的差分干涉图进行人工标注,建立神经网络模型的训练集和测试集。通过旋转、偏移、缩放、镜像等数据增强方法扩大数据集,最后以4:1的比例划分训练集和测试集。模型训练时,通过观察损失函数的收敛情况调整神经网络模型的学习率、进入模型的批样本数量(batch-size)、动量参数和权重衰减系数,使损失函数达到收敛。最后,对S3处理的所有差分干涉图进行预测,获取全部的采空区识别信息。
S5、对识别后采空区-非采空区二值图像栅格转矢量,生成矢量化的采空区结果。
对识别后采空区图像栅格转矢量,具体步骤为:将U-net提取的二分类栅格图像(即采空区和非采空区),根据二分类图像的梯度变化,将栅格采空区图像进行矢量化处理,生成矢量化的采空区结果。
经过S4预测后的图像是一个二分类的栅格图像,只存在采空区和非采空区两种信息。为了进行后续的矢量融合和时空分析,利用二分类栅格图像梯度变化信息,对栅格图像矢量化,提取采空区的矢量信息。
S6、对多期矢量化后的采空区进行空间融合。
所述的对多期矢量化后的采空区结果进行空间融合,具体步骤为:针对多期矢量化后的采空区结果,用空间拓扑分析进行矢量融合,将每期的小范围矢量取并集融合成整体采空区矢量。
以短时空基线的准则建立差分干涉对后,得到多期不同时间跨度的差分干涉对,每期干涉对均提取了采空区的矢量信息,其采空区的形状、大小、分布特点均有不同。为了得到最终的采空区边界范围,需要对多期采空区矢量进行范围取并集的融合(见图1)。在融合前先对各期采空区矢量建立时间、面积、质心属性字段信息,融合后将所有字段保留,便于后续的时空分析。
S7、时空分析采空区面积大小、发展方向和采空区边界距铁路的最近距离。
统计分析采空区面积大小、发展方向和采空区边界距铁路的最近距离,具体步骤为:对融合后的整体采空区矢量进行统计分析,获取采空区面积大小,发展方向以及采空区边缘距铁路的最近距离。
经过S6获取的采空区矢量是最终的整体采空区矢量边界,利用ArcGIS添加面积字段属性,统计分析每个采空区的面积大小。
针对每期识别的采空区矢量,利用ArcGIS计算其矢量的质心坐标,通过多项式函数,以最小二乘为准则,拟合质心变化曲线。曲线拟合公式如下所示:
p(x)=a0+a1x+a2x2+…an*xk
Figure BDA0002497724850000051
式中:xi为质心x坐标,yi为质心y坐标
每个质心点带有时间属性信息,按照时间出现的先后顺序作为曲线的增长方向,以此作为开采发展方向。
同时,借助外部铁路线位信息,在ArcGIS中将线位与开采区矢量图层进行叠加,生成近邻表得到每个采空区矢量与线位的最邻近距离,并将最邻近距离添加到矢量对应的字段属性信息。
最后将所有影响因子和外部铁路线位矢量文件叠加到ArcGIS中,系统分析采空区分布特点、发展方向、面积范围、边缘距铁路最近距离信息,分析评价采空区对铁路的综合影响,形成最终报告。
显然,上述实施例仅是为清楚地说明所作的举例,而并非对实施方式的限定。对于所属领域的普通技术人员来说,在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。

Claims (9)

1.一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于,包括以下步骤:
(1)获取研究区域内时间序列SAR影像和精密轨道文件;
(2)将SAR影像数据预处理建立干涉对连接图;
(3)对连接图中的干涉对进行干涉处理生成差分干涉图,将差分干涉图进行地理编码;
(4)利用U-net神经网络对地理编码后差分干涉图进行采空区识别生成采空区-非采空区二值图像;
(5)对识别后采空区-非采空区二值图像由栅格转矢量,生成矢量化的采空区结果;
(6)对多期矢量化后的采空区结果进行空间融合;
(7)统计分析采空区面积大小、发展方向和采空区边界距铁路的最近距离。
2.如权利要求1所述的一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于:获取研究区域内时间序列SAR影像和精密轨道文件,具体步骤为:
每月时间段内所获取的同极化SLC影像不少于2幅;同时,获取对应日期SAR影像的卫星精密轨道参数文件。
3.如权利要求1所述的一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于:将SAR影像数据预处理建立连接图,具体步骤为:
将所有SAR影像聚焦、配准到超级主影像,对影像进行去斜,用精密轨道文件优化影像参数文件中的轨道参数,根据矢量文件裁剪覆盖研究区的范围;再对所有影像建立干涉对连接图,将未自动参与组网影像按时空最邻近原则手动添加参与组网。
4.如权利要求3所述的一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于:影像空间垂直基线为临界基线的0-30%,时间基线为0-95d,每幅影像干涉对数应至少保证2对。
5.如权利要求1所述的一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于:所述的干涉处理具体步骤为:
对建立连接图的影像进行差分干涉、去平、轨道精炼、重去平、滤波处理;所述的地理编码差分干涉图是将差分干涉图从雷达坐标系转换到WGS84地理坐标系。
6.如权利要求1所述的一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于:利用U-net神经网络对地理编码后图像进行采空区识别生成采空区-非采空区二值图像,具体步骤为:
利用U-net神经网络结构,对编码后的差分干涉图进行语义分割,对采空区图像进行像素级别的分类,从而进行采空区识别生成采空区-非采空区二值图像。
7.如权利要求1所述的一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于:对识别后采空区-非采空区二值图像栅格转矢量,具体步骤为:
对U-net提取的栅格二值图像(即采空区和非采空区),根据二分类图像的梯度变化,将栅格采空区图像进行矢量化处理,生成矢量化的采空区结果。
8.如权利要求1所述的一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于:所述的对多期矢量化后的采空区结果进行空间融合,具体步骤为:
针对多期矢量化后的采空区结果,用空间拓扑分析进行矢量融合,将每期的小范围矢量取并集融合成整体采空区矢量。
9.如权利要求1所述的一种基于D-InSAR的铁路采空区边界提取新方法,其特征在于:统计分析采空区面积大小、发展方向和采空区边界距铁路的最近距离,具体步骤为:
对融合后的整体采空区矢量进行统计分析,获取采空区面积大小,发展方向以及采空区边缘距铁路的最近距离。
CN202010423214.5A 2020-05-19 2020-05-19 一种基于D-InSAR的铁路采空区边界提取新方法 Active CN111623749B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010423214.5A CN111623749B (zh) 2020-05-19 2020-05-19 一种基于D-InSAR的铁路采空区边界提取新方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010423214.5A CN111623749B (zh) 2020-05-19 2020-05-19 一种基于D-InSAR的铁路采空区边界提取新方法

Publications (2)

Publication Number Publication Date
CN111623749A CN111623749A (zh) 2020-09-04
CN111623749B true CN111623749B (zh) 2021-11-16

Family

ID=72270051

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010423214.5A Active CN111623749B (zh) 2020-05-19 2020-05-19 一种基于D-InSAR的铁路采空区边界提取新方法

Country Status (1)

Country Link
CN (1) CN111623749B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112184902B (zh) * 2020-09-21 2022-09-02 东华理工大学 一种面向越界开采识别的地下开采面反演方法
CN114063063A (zh) * 2021-09-09 2022-02-18 浙江浙能天然气运行有限公司 基于合成孔径雷达和点状传感器的地质灾害监测方法
CN115331096B (zh) * 2022-06-28 2023-04-21 中国自然资源航空物探遥感中心 一种采矿沉陷区域识别方法、系统、存储介质和电子设备

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106226764A (zh) * 2016-07-29 2016-12-14 安徽理工大学 一种基于D‑InSAR的煤矿开采地沉陷区域的测定方法
IL258119B2 (en) * 2018-03-14 2024-06-01 Elta Systems Ltd Coherent change detection techniques

Also Published As

Publication number Publication date
CN111623749A (zh) 2020-09-04

Similar Documents

Publication Publication Date Title
CN111623749B (zh) 一种基于D-InSAR的铁路采空区边界提取新方法
Papasaika et al. Fusion of digital elevation models using sparse representations
CN109100719B (zh) 基于星载sar影像与光学影像的地形图联合测图方法
CN104574512A (zh) 一种顾及地形语义信息的多尺度dem构建方法
CN109489625A (zh) 一种城市区域地表形变监测方法
CN109212522A (zh) 一种获取数字地图的方法和装置
CN114325706B (zh) 一种分布式散射体滤波方法
CN114966689A (zh) 沿海城市时序InSAR沉降监测分析方法、装置、设备及介质
CN113251914A (zh) InSAR技术与长短时记忆神经网络结合的地表形变预测方法
CN114791273A (zh) 一种针对滑坡的InSAR形变监测结果解释方法
Rastogi et al. Urban sprawl analysis using shannon’s entropy and fractal analysis: A case study on Tiruchirappalli city, India
Cheng et al. Urban land extraction using DMSP/OLS nighttime light data and OpenStreetMap datasets for cities in China at different development levels
CN118196637A (zh) 一种基于InSAR形变与影响因子耦合的潜在滑坡识别方法
CN112097733A (zh) InSAR技术与地理探测器结合的地表形变监测方法
CN115638767A (zh) 地面沉降监测方法
Luo et al. Dynamic analysis of urban ground subsidence in Beijing based on the permanent scattering InSAR technology
CN113238228B (zh) 基于水准约束的三维地表形变获取方法、系统及装置
CN117236518B (zh) 一种地铁沿线地质形变的预测方法
CN114170524A (zh) 一种高时-空卫星数据融合的单日潮滩数字地形构建方法
CN116863029A (zh) 一种基于多源遥感数据的地震灾害应急制图方法及装置
Yang et al. Ground deformation monitoring of Zhengzhou city from 2012 to 2013 using an improved IPTA
Hapep et al. Comparison of Different DEM Generation Methods based on Open Source Datasets.
Amhar Quality test various existing DEM in Indonesia toward 10 meter national DEM
CN114114256A (zh) 一种基于D-InSAR-GIS叠加分析技术的大面积矿区沉陷监测方法
Khemiri et al. DEM-Based GIS Algorithms and 3D Geospatial Mapping for Creation of Hydrogeological Models Data in Foussana Basin (Central Tunisia)

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