CN114417724B - 一种山地城市土地利用演变的模拟方法 - Google Patents

一种山地城市土地利用演变的模拟方法 Download PDF

Info

Publication number
CN114417724B
CN114417724B CN202210069578.7A CN202210069578A CN114417724B CN 114417724 B CN114417724 B CN 114417724B CN 202210069578 A CN202210069578 A CN 202210069578A CN 114417724 B CN114417724 B CN 114417724B
Authority
CN
China
Prior art keywords
land utilization
simulation
land
mountain
conversion
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
CN202210069578.7A
Other languages
English (en)
Other versions
CN114417724A (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
Original Assignee
Hohai University HHU
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 filed Critical Hohai University HHU
Priority to CN202210069578.7A priority Critical patent/CN114417724B/zh
Publication of CN114417724A publication Critical patent/CN114417724A/zh
Application granted granted Critical
Publication of CN114417724B publication Critical patent/CN114417724B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/243—Classification techniques relating to the number of classes
    • G06F18/24323—Tree-organised classifiers
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种山地城市土地利用演变的模拟方法，包括如下具体步骤：S1、预处理数据：获取山地城市遥感影像数据，并进行预处理，S2、确定驱动因子，贡献分析：选取该城市驱动因子数据，并对不同土地利用类型进行因子贡献分析，S3、模型模拟：通过确定转换规则、邻域权重等参数，对土地利用情况进行模拟，S4、精度评价：将模拟结果与真实土地利用数据对比，通过Kappa系数和FoM指数验证精度。本发明通过在改进后的CA模型‑PLUS模型的基础上进行模拟，该模型中新的转换规则数据挖掘框架可以用于识别土地利用变化的规则，结合了转化分析策略和格局分析策略的优势，通过增加山地城市特殊的驱动因子，区分各因子贡献度，提高山地城市的土地利用模拟精度，从而精准分析山地城市土地利用演变的驱动因子，精确模拟土地利用变化。

Description

一种山地城市土地利用演变的模拟方法
技术领域
本发明涉及地理信息技术领域，具体为一种山地城市土地利用演变的模拟方法。
背景技术
目前，现有的土地利用模拟模型采用较多的是元胞自动机模型及衍生的一系列改进模型，为了提高模拟的精度和运算速度，改进的元胞自动机模型技术逐渐发展，如在申请号为201410409993.8的中国专利“一种基于随机森林的元胞自动机城市扩展模拟方法”中，提出随机森林-元胞自动机模型；
但是目前对城市土地利用演变的模拟中，对于山地等特殊城市的模拟较少，山地城市由于其地势较高、落差较大、坡度大，各土地利用类型的分布会存在明显受地势影响的情况，且由于地势、山地阻隔原因，城镇、组团间相互连通的道路网络等建设用地的布设也深受影响，因此，对于山地城市土地利用的模拟分析显得尤为重要。
发明内容
本发明提供一种山地城市土地利用演变的模拟方法，可以有效解决上述背景技术中提出目前对城市土地利用演变的模拟中，对于山地等特殊城市的模拟较少，山地城市由于其地势较高、落差较大、坡度大，各土地利用类型的分布会存在明显受地势影响的情况，且由于地势、山地阻隔原因，城镇、组团间相互连通的道路网络等建设用地的布设也深受影响的问题。
为实现上述目的，本发明提供如下技术方案：一种山地城市土地利用演变的模拟方法，包括如下具体步骤：
S1、预处理数据：获取山地城市遥感影像数据，并进行预处理；
S2、确定驱动因子，贡献分析：选取该城市驱动因子数据，并对不同土地利用类型进行因子贡献分析；
S3、模型模拟：通过确定转换规则、邻域权重等参数，对土地利用情况进行模拟；
S4、精度评价：将模拟结果和真实土地利用数据进行对比，通过Kappa系数和FoM指数验证精度。
根据上述技术方案，所述S2中，分析各因子对于不同土地利用类型变化的贡献，采用均方根误差和随机森林特有的袋外均方根误差检测值计算自变量对因变量变化的重要性；
采用随机森林方法对各个地类进行随机采样，可以用于输出i栅格的土地利用类型k的增长概率公式如下：
式(1)中，d＝1表示存在由其他用地类型转换成k类型，其余情况用d＝0表示；
x表示驱动因子向量；
I(·)函数是关于随机森林中的决策树集；
hn(x)表示向量x在第n个决策树下的预测类型；
m是所有决策树的个数；
结合各土地利用类型的发展概率、邻域权重、转换规则以及区域约束条件实现对未来各土地利用类型像元空间分布的合理化配置。
根据上述技术方案，所述S3中，利用PLUS模型进行模拟，模拟时元胞数量变化，在模拟时需对模型中的邻域权重参数反复的进行调试，直至模拟精度达到要求。
根据上述技术方案，所述S4中，Kappa系数综合了用户精度和制图精度两个参数来检验栅格土地利用类型空间精度，计算公式为：
式(2)中，Ps为实际土地利用类型与模拟的土地利用类型相一致的比例，即模拟的正确率；
Pr为随机条件下期望的正确率；
Kappa系数的取值范围是在0～1之间，Kappa系数越接近于1，说明一致性越好，模拟的土地利用图与实际土地利用图越接近。
FoM指数衡量实际观测转换数量与模拟预测转换之间的一致性，具体计算公式为：
式(3)中，NA表示实际观测中发生转换而模拟过程中未发生转换的像元数；
NB表示实际观测中发生转换且模拟正确的像元数；
NC表示实际观测中发生变化而模拟过程中发生转换但转换类型错误的像元数；
ND表示实际观测中未发生转换而模拟过程中发生转换的像元数。
根据上述技术方案，所述S1中，获取起始和终止年份遥感影像并进行预处理，结合当地山地城市土地利用/覆被实际情况，进行土地利用分类。
根据上述技术方案，所述S2中，筛选山地城市驱动因子原始数据，除一般GDP、人口密度驱动因子外，增加高程、坡度以及多级道路网路数据，并进行欧氏距离测算。
与现有技术相比，本发明的有益效果：
1、通过在改进后的CA模型-PLUS模型的基础上进行模拟，该模型中新的转换规则数据挖掘框架可以用于识别土地利用变化的规则，结合了转化分析策略和格局分析策略的优势，通过增加山地城市特殊的驱动因子，区分各因子贡献度，提高山地城市的土地利用模拟精度，从而精准分析山地城市土地利用演变的驱动因子，精确模拟土地利用变化。
2、考虑山地城市高程和道路的特点，采用改进的元胞自动机PLUS模型，基于土地扩张分析策略的转化规则挖掘，在模拟中增加多级道路网络、高程和坡度因子，提高模拟精度，较好的分析识别了山地城市高程、坡度以及多级道路网路等各驱动因子的贡献，为后续技术提升、科学研究奠定了基础，同时为国家对于山地城市土地利用布局与发展提供了实际决策依据。
附图说明
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。
在附图中：
图1是本发明山地城市土地利用演变的模拟步骤图；
图2是本发明的技术流程图；
图3是本发明各土地利用类型发展概率图；
图4是本发明各土地利用类型的用地扩张因子贡献图；
图5是本发明各土地利用类型的训练精度；
图6是本发明模拟时元胞数量变化图；
图7是本发明两江新区2015年土地利用模拟结果。
具体实施方式
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
实施例1：
如图1-2所示，本发明提供一种技术方案，一种山地城市土地利用演变的模拟方法，包括如下具体步骤：
S1、预处理数据：获取山地城市遥感影像数据，并进行预处理；
S2、确定驱动因子，贡献分析：选取该城市不同类型的地块和驱动因子，并进行贡献分析；
S3、模型模拟：在模型中输入参数，调试参数对土地的利用情况进行模拟；
S4、精度评价：将模拟结果和真实土地利用数据进行对比，通过Kappa系数和FoM指数测算模拟精度。
根据上述技术方案，S2中，分析各因子对于不同土地利用类型变化的贡献，采用均方根误差和随机森林特有的袋外均方根误差检测值计算自变量对因变量变化的重要性；
采用随机森林方法对各个地类进行随机采样，可以用于输出i栅格的土地利用类型k的增长概率公式如下：
式(1)中，d＝1表示存在由其他用地类型转换成k类型，其余情况用d＝0表示；
x表示驱动因子向量；
I(·)函数是关于随机森林中的决策树集；
hn(x)表示向量x在第n个决策树下的预测类型；
m是所有决策树的个数；
结合各土地利用类型的发展概率、邻域权重、转换规则以及区域约束条件实现对未来各土地利用类型像元空间分布的合理化配置。
根据上述技术方案，S3中，利用PLUS模型进行模拟，模拟时元胞数量变化，在模拟时需对模型中的邻域权重参数反复的进行调试，直至模拟精度达到要求。
根据上述技术方案，S4中，Kappa系数综合了用户精度和制图精度两个参数来检验栅格土地利用类型空间精度，计算公式为：
式(2)中，Ps为实际土地利用类型与模拟的土地利用类型相一致的比例，即模拟的正确率；
Pr为随机条件下期望的正确率；
Kappa系数的取值范围是在0～1之间，Kappa系数越接近于1，说明一致性越好，模拟的土地利用图与实际土地利用图越接近。
FoM指数衡量实际观测转换数量与模拟预测转换之间的一致性，具体计算公式为：
式(3)中，NA表示实际观测中发生转换而模拟过程中未发生转换的像元数；
NB表示实际观测中发生转换且模拟正确的像元数；
NC表示实际观测中发生变化而模拟过程中发生转换但转换类型错误的像元数；
ND表示实际观测中未发生转换而模拟过程中发生转换的像元数。
根据上述技术方案，S1中，获取起始和终止年份遥感影像并进行预处理，结合当地山地城市土地利用/覆被实际情况，进行土地利用分类。
根据上述技术方案，S2中，筛选山地城市驱动因子原始数据，除一般GDP、人口密度驱动因子外，增加高程、坡度以及多级道路网路数据，并进行欧氏距离测算。
实施例2：
以重庆市两江新区为例，实施步骤如下：
步骤一：预处理数据
在Google Earth Engine平台获取两期两江新区遥感影像，结合Google Earth高清影像，通过目视解释对影像进行分类，获取起始年份和终止年份土地利用数据，分类体系结合当地山地城市具体情况，分为：建设用地(主要包括城镇工矿用地、城镇居民用地、交通水利用地)、植被(主要包括乔木、灌木多种森林植被)、裸地(主要包括未利用地)和水体(主要包括河流、湖泊水库)。
在OpenStreetMap平台、中科院资源环境科学与数据中心、重庆市地理信息中心、地理空间数据云平台获取基础驱动因子数据，使用ArcGIS计算两江新区计算栅格点到各要素的空间距离，得到表1数据。
表1两江新区驱动因子选择及其含义
步骤二：确定驱动因子，贡献分析
如图3所示，利用随机采样的方式对植被、建设用地、水体和裸地四类土地利用类型进行采样，采样率设置为0.01，决策树的数量设为12，最后得到两江新区2015年土地利用类型在各像元上的发展概率和各驱动因子的贡献度。
如图4所示，从各土地利用类型用地扩张因子贡献结果来看，影响两江新区植被分布的因素主要有到火车站的距离、到二级道路的距离、GDP且与当地的GDP呈明显的正相关，植被分布同与火车站的距离以及到二级道路的距离呈显著的负相关，重要性均大于0.11，结果表明，交通便利、经济发达的地方植被分布较多，反之亦然，影响两江新区建设用地分布的因素主要有高程、到二级道路的距离、到火车站的距离要素，且与高程、到二级道路的距离、到火车站的距离均为负相关。表明建设用地多分布于交通便利的平坦地区，高程较高、坡度较大的地方不利于房屋建设，且建设用地通常考虑到居住的舒适性，所以与距道路一定距离修建的事实相符，而二级道路对当地居民的生产生活影响较大。水体分布主要与高程呈负相关，高程驱动力因子的重要性高达0.502，表明两江新区的水体，主要分布在地势低平的地区。影响裸地分布的因素主要有到二级道路的距离、到铁路的距离、到一级道路的距离要素，且均呈负相关。
如图5所示，采用最常见的均方根误差和随机森林特有的袋外均方根误差检测值对预测结果进行精度评定，训练结果中，水体是四种土地利用中误差最小的，经检验，两江新区2015年土地利用训练的精度较高且误差在允许的范围之内。
步骤三：模型模拟
在模拟时需对模型中的邻域权重参数反复的进行调试，直至模拟精度达到要求，最终，各权重因子设置为：植被0.35、建设用地0.95、水体0.4、裸地0.2；
考虑到两江新区内的长江和嘉陵江，是重庆最重要的河流水域，它们不仅是当地及下游居民的水源保证，也是水运交通要道，因此对研究区域内水体类型进行限制，用二值化的栅格图来表示，其中不允许水体发生变化的区域设值为0，允许水体发生变化的区域设置为1；
如图6所示，利用PLUS模型进行模拟，模拟时元胞数量变化；
如图7所示，2015年土地利用模拟结果。
步骤四：精度评价
将2015年模拟图与实际土地利用图利用栅格计算器进行相比，得到模拟正确的值为1，统计发现模拟正确的像元数为1077694个，实际总像元数为1310414个，则模拟精度为82％，理想状态下正确模拟率1，根据公式计算得到2015年两江新区模拟图的Kappa指数为0.73，FoM指数为0.263，说明此次模拟效果是具有可行性的。
最后应说明的是：以上所述仅为本发明的优选实例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种山地城市土地利用演变的模拟方法，其特征在于：包括如下具体步骤：
S1、预处理数据：获取山地城市遥感影像数据，并进行预处理；
S2、确定驱动因子，贡献分析：选取该城市驱动因子数据，并对不同土地利用类型进行因子贡献分析；
S3、模型模拟：通过确定转换规则、邻域权重参数，对土地利用情况进行模拟；
S4、精度评价：将模拟结果和真实土地利用数据进行对比，通过Kappa系数和FoM指数验证精度；
所述S2中，分析各因子对于不同土地利用类型变化的贡献，采用均方根误差和随机森林特有的袋外均方根误差检测值计算自变量对因变量变化的重要性；
采用随机森林方法对各个地类进行随机采样，用于输出i栅格的土地利用类型k的增长概率公式如下：
式(1)中，d＝1表示存在由其他用地类型转换成k类型，其余情况用d＝0表示；
x表示驱动因子向量；
I(·)函数是关于随机森林中的决策树集；
hn(x)表示向量x在第n个决策树下的预测类型；
m是所有决策树的个数；
结合各土地利用类型的发展概率、邻域权重、转换规则以及区域约束条件实现对未来各土地利用类型像元空间分布的合理化配置；
所述S4中，Kappa系数综合了用户精度和制图精度两个参数来检验栅格土地利用类型空间精度，计算公式为：
式(2)中，Ps为实际土地利用类型与模拟的土地利用类型相一致的比例，即模拟的正确率；
Pr为随机条件下期望的正确率；
Kappa系数的取值范围是在0～1之间，Kappa系数越接近于1，说明一致性越好，模拟的土地利用图与实际土地利用图越接近；
FoM指数衡量实际观测转换数量与模拟预测转换之间的一致性，具体计算公式为：
式(3)中，NA表示实际观测中发生转换而模拟过程中未发生转换的像元数；
NB表示实际观测中发生转换且模拟正确的像元数；
NC表示实际观测中发生变化而模拟过程中发生转换但转换类型错误的像元数；
ND表示实际观测中未发生转换而模拟过程中发生转换的像元数。
2.根据权利要求1所述的一种山地城市土地利用演变的模拟方法，其特征在于，所述S3中，利用PLUS模型进行模拟，模拟时元胞数量变化，在模拟时需对模型中的邻域权重参数反复的进行调试，直至模拟精度达到要求。
3.根据权利要求1所述的一种山地城市土地利用演变的模拟方法，其特征在于，所述S1中，获取起始和终止年份遥感影像并进行预处理，结合当地山地城市土地利用/覆被实际情况，进行土地利用分类。
CN202210069578.7A 2022-01-20 2022-01-20 一种山地城市土地利用演变的模拟方法 Active CN114417724B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210069578.7A CN114417724B (zh) 2022-01-20 2022-01-20 一种山地城市土地利用演变的模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210069578.7A CN114417724B (zh) 2022-01-20 2022-01-20 一种山地城市土地利用演变的模拟方法

Publications (2)

Publication Number Publication Date
CN114417724A CN114417724A (zh) 2022-04-29
CN114417724B true CN114417724B (zh) 2024-07-12

Family

ID=81275555

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210069578.7A Active CN114417724B (zh) 2022-01-20 2022-01-20 一种山地城市土地利用演变的模拟方法

Country Status (1)

Country Link
CN (1) CN114417724B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114861277B (zh) * 2022-05-23 2022-11-15 中国科学院地理科学与资源研究所 一种长时序国土空间功能与结构模拟方法
CN115600075B (zh) * 2022-12-12 2023-04-28 深圳市城市规划设计研究院股份有限公司 景观斑块变化的测度方法和装置、电子设备及存储介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111080070A (zh) * 2019-11-19 2020-04-28 同济大学 一种基于空间误差的城市土地利用模拟元胞自动机方法
CN112131731A (zh) * 2020-09-15 2020-12-25 同济大学 一种基于空间特征向量滤波的城市生长元胞模拟方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106295576B (zh) * 2016-08-12 2017-12-12 中国水利水电科学研究院 一种基于自然地理特征的水源类型解析方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111080070A (zh) * 2019-11-19 2020-04-28 同济大学 一种基于空间误差的城市土地利用模拟元胞自动机方法
CN112131731A (zh) * 2020-09-15 2020-12-25 同济大学 一种基于空间特征向量滤波的城市生长元胞模拟方法

Also Published As

Publication number Publication date
CN114417724A (zh) 2022-04-29

Similar Documents

Publication Publication Date Title
CN110472882B (zh) 基于主成分分析的城市开发用地适宜性评价方法
CN114417724B (zh) 一种山地城市土地利用演变的模拟方法
CN112163367B (zh) 融合萤火虫算法和元胞自动机的城市扩张模拟预测方法
CN111598314B (zh) 一种基于绩效单元的公园绿地规划布局方法
CN112884226A (zh) 基于多智能体算法的国土空间格局模拟规划方法及***
CN110909924A (zh) 基于交叉熵优化器的城市扩张多情景模拟元胞自动机方法
CN110889196B (zh) 基于水质模型的水环境承载力评估方法、装置及存储介质
US20150204922A1 (en) Method for Predicting Wind Power Density
CN110826244B (zh) 模拟轨道交通影响城市生长的共轭梯度元胞自动机方法
CN106845559A (zh) 顾及poi数据空间异质性的地表覆盖验证方法及***
CN110134907B (zh) 一种降雨缺失数据填补方法、***及电子设备
CN113191582B (zh) 一种基于gis与机器学习的道路山洪易损性评价方法
CN112784503B (zh) 一种闸坝调度下水鸟栖息地适宜性计算方法与装置
CN115271373A (zh) 城市群弹性开发边界划定方法及***
CN114398951A (zh) 一种基于随机森林和众源地理信息的土地利用变化驱动因子挖掘方法
CN113298439A (zh) 基于人口分布的环境风险评估方法、装置及计算机设备
CN116415110A (zh) 基于多源遥感密度数据进行碳排放量分区网格化的方法
CN113011720B (zh) 一种基于ca模型的增量递减的城镇开发边界划定方法
CN114219370A (zh) 一种基于社交网络的河流水质多维影响因素权重分析方法
Momeni et al. Pattern‐based calibration of cellular automata by genetic algorithm and Shannon relative entropy
Badhrudeen et al. Machine learning based tool for identifying errors in CAD to GIS converted data
CN113743659A (zh) 一种基于成分法和马尔可夫元胞自动机的城市布局预测方法及应用
CN117709580A (zh) 一种基于setr与地理网格的海洋承灾体脆弱性评价方法
CN115829163B (zh) 基于多模式集成的长江中下游流域径流预测方法和***
CN115619275B (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