CN110310370B - 一种gps与srtm点面融合的方法 - Google Patents

一种gps与srtm点面融合的方法 Download PDF

Info

Publication number
CN110310370B
CN110310370B CN201910528151.7A CN201910528151A CN110310370B CN 110310370 B CN110310370 B CN 110310370B CN 201910528151 A CN201910528151 A CN 201910528151A CN 110310370 B CN110310370 B CN 110310370B
Authority
CN
China
Prior art keywords
grid
value
elevation
point
gps
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
CN201910528151.7A
Other languages
English (en)
Other versions
CN110310370A (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.)
Central South University of Forestry and Technology
Original Assignee
Central South University of Forestry and Technology
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 Central South University of Forestry and Technology filed Critical Central South University of Forestry and Technology
Priority to CN201910528151.7A priority Critical patent/CN110310370B/zh
Publication of CN110310370A publication Critical patent/CN110310370A/zh
Application granted granted Critical
Publication of CN110310370B publication Critical patent/CN110310370B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Graphics (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

一种基于积分平差实现GPS与SRTM点面融合的方法，包括以下步骤：1)通过地面GPS观测，获得目标地区SRTM的DEM格网上任一点p的GPS高程值hp；2)用测量平差方法把hp和p所在DEM格网的高程值
Figure DDA0002098861800000011
进行融合；得到p点所在格网的格网值和四个角点a、b、c、d新的高程值；作为第一层融合；3)以a、b、c、d四个点为新的观测值，分别对其周围的网格进行融合，得到外围8格网融合后新的格网值及外围12个节点的新的高程；作为第二层融合；依次类推，得到全格网的高程值；当有格网的角点得到多个高程值时，取平均值。本发明的基于积分平差实现GPS与SRTM点面融合的方法的精度高。

Description

一种GPS与SRTM点面融合的方法
技术领域
本发明涉及一种数字高程模型的处理方法，尤其涉及一种基于积分平差实现GPS与 SRTM点面融合的方法。
背景技术
数字高程模型(Digital Elevation Model，简称DEM)是通过有限的地形高程数据实现对地形曲面的数字化模拟，它是对二维地理空间上具有连续变化特征地理现象的模型化表达和过程模拟。如今，DEM已被广泛应用，如土木工程、军事应用、地质灾害和可视化分析等。
建立DEM实现地表重建的前提是对连续地球表面进行有限离散点的精确采样和三维坐标测量。目前，DEM数据获取方法主要包括野外测量、摄影测量和遥感、制图数字化、合成孔径雷达和机载激光扫描。如今，使用卫星图像生成DEM有一个显著的优势因为它相对便宜而且生成DEM所需的时间更少。然而使用光谱范围的缺点是它需要高分辨率，良好的光照条件和无云视线才能得到较好的DEM精度。
近年来，干涉合成孔径雷达(InSAR)因其不依赖于自然光照而作为一个主动系统在提取高程数据方面变成流行技术。2000年美国航天飞机雷达地形测绘任务(Shuttleradar topography mission,SRTM)利用机载Insar技术，仅11d便获取了全球80％陆地1″分辨率的三维地形信息。2014年9月，1″分辨率的SRTM数据逐步向全球用户免费开放。SRTM 在平坦区域的准确度范围为2.9-3.4米，半平坦区域为5.4-6.29米([1]A.Karwel andI.Ewiak, "Estimation of the accuracy of the SRTM terrain model on the area ofPoland,"The international archives of the photogrammetry,remote sensing andspatial information sciences,vol.37,no.B7,pp.169-172,2008.)。
已经有很多方法用来提高SRTM的精度，例如融合Aster和SRTM数据、融合ICESat数据和SRTM数据等；在我国，每个省市都建有CORS站，每天收集全省各地大量点的GPS 高程(或北斗高程)，充分利用这些GPS高程数据具有重要的实际和实用价值。
发明内容
本发明要解决的技术问题是克服现有技术的不足，提供一种精度高的基于积分平差实现 GPS与SRTM点面融合的方法。
为解决上述技术问题，本发明提出的技术方案为：
如图1所示，某地区SRTM的DEM，假设任意一点p所在DEM格网的高程值为
Figure BDA0002098861780000011
通过地面GPS观测，获得p点的GPS高程值hp，我们的目标是把p点的GPS高程值hp与该地区的SRTMDEM进行融合，以达到提高DEM的精度。融合的基本思想是，首先用测量平差的方法把hp与
Figure BDA0002098861780000023
进行融合，求得p点所在格网的格网值和四个角点a,b,c,d新的高程值。此为以p点为中心的中心融合。a,b,c,d四个点的高程又会分别影响与它们相邻外围的8个相邻网格的高程，因此可以把这4个点的高程看作它们外围相邻格网的新观测值，分别与相邻格网进一步融合，求得外围8个格网融合后新的格网值及外围12个节点的新的高程。此为以p点为中心的第二层融合。外围12个节点的新的高程又会影响围绕第二层的16个格网，逐个格网融合，就可完成以p点为中心的第三层融合。依此类推。p点的高程值hp通过这个方式可以提高其他格网的高程值精度。由于一个格网会有多个相邻格网，有些格网的角点会得到多个新值，可取平均作为最后值。
以图1为例，在第二层融合过程中，格网(i-1,j)与格网(i-1,j+1)有共同的节点，里面的节点a的高程使用第一层的结果，外围节点的高程值则取第二层融合时，相邻两格网计算结果的平均值。
具体的平差模型可按如下方法建立：
一般而言，DEM中每个格网的高程值
Figure BDA0002098861780000024
表示这个格网内所有点高程的平均值。用h(x,y) 表示任意一点(x,y)的高程，
Figure BDA0002098861780000025
表示格网(i,j)的格网值,那么
Figure BDA0002098861780000026
可以用h(x,y)在这个格网内的积分的平均值表示，即：
Figure BDA0002098861780000021
其中，S为格网面积。DEM的各格网值一般是通过各种测绘手段获取的
Figure BDA0002098861780000027
的观测值hi,j。考虑观测误差，有：
Figure BDA0002098861780000022
p点的高程可表示为h(xp，yp)，考虑到观测误差的存在，p点高程的观测值可表示为：
hp+v＝h(xp，yp) (3)
因为实际地形一般都非常复杂，故难以用数学函数准确描述，所以很难求得h(x，y)的精确表述。我们可以考虑一个非常小的局部，比方说，一个围绕格网(i,j)的3*3的格网范围，在这个小区域内，找一个数学函数来近似代替h(x，y)是可行的，例如可用一个二次函数代替h(x，y)等。利用这个区域9个格网点的观测值与P点的观测值，用测量平差的方法可联合求解这个区域的数学函数(一个高程面估计)，根据这个高程面，可以求得融合后格网(i，j)和四个角点的高程值，再由式(1)可以重新确定融合后格网(i，j)的高程。按式(1)可构成这个区域9个点的观测方程，再与式(2)联合，可形成如下平差模型：
Figure BDA0002098861780000031
hp+vp＝h(x，y) (5)
如果一个二次函数代替h(x，y)，就有
h(x，y)＝a1+a2x+a3y+a4xy+a5x2+a6y2 (6)
Figure BDA0002098861780000032
hp+v＝a1+a2xp+a3yp+a4xpyp+a5xp 2+a6yp 2 (8)
最后可构成平差模型：
L+V＝AX (9)
其中
L＝(hi-1，j-1 hi-1，j...hi+1，j hi+1，j+1 hp)T
V＝(Vi-1，j-1 Vi-1，j...Vi+1，j Vi+1，j+1 Vp)T
X＝(a1 a2 a3 a4 a5 a6)T
Figure BDA0002098861780000033
根据最小二乘平差，有
Figure BDA0002098861780000035
其中P为权值，定权为900。平差融合后，小区域内各点的高程值可表示为
Figure BDA0002098861780000034
格网(i,j)平差融合后的值分别为：
Figure BDA0002098861780000041
格网(i,j)四个角点的融合后的高程分别为
Figure BDA0002098861780000042
Figure BDA0002098861780000043
Figure BDA0002098861780000044
Figure BDA0002098861780000045
格网(i,j)平差融合后的值
Figure BDA0002098861780000046
即为该格网融合后的格网值。格网(i,j)四个角点融合后的高程值
Figure BDA0002098861780000047
对于相邻的格网就是新的观测值，如对于格网(i-1,j)，
Figure BDA0002098861780000048
是2个新的观测，新观测值的权定为900。按与上述一样的方法，取格网(i-1,j)及周边9个格网值与这2个新观测一起进行平差，可以求得该格网融合后的高程值及格网4个角点融合后的高程值。具体可按如下步骤进行：
(1)取p点及所在格网与周边8个格网点一起平差，定权为900，求得所在格网及4个角点融合后的高程值，融合后的格网值作为该格网最后的高程值，此为以p点为中心的第一层融合。
(2)经过第一层融合后，围绕格网(i,j)的8个网格会受到格网(i,j)的4个角点高程变化的影响。其中格网(i-1,j)有a、c两个新观测点，格网(i,j+1)有a、b两个新观测点，格网(i+1,j)有b、c两个新观测，格网(i,j-1)有c、d两个新观测点。格网(i-1,j+1) 有新观测点a，格网(i+1,j+1)有新观测点b，格网(i+1,j)有新观测点d，格网(i-1,j-1) 有新观测点c。按上述平差方法分别将新观测点的观测与相应格网融合，可求得第二层8个格网的格网融合后的新值，及外围节点的新的高程值。对于相邻网格共同的节点取两个网格各自计算出的节点新高程值的平均值作为最后的新值。
(3)计算第二层8个网格平差融合后的新值与原格网值的差值，若各差值的均值<阀值 (0.5米)，则认为整个改进效果不明显，计算终止，算法收敛。对于各差值的均值>阀值(0.5 米)的格网，表明原DEM有误差需要修正，则与该格网相关的计算可以进一步往下走进行，直至算法收敛。
本发明提出的实现GPS与SRTM点面融合的积分平差算法，该算法的特点是把任意一点的高程表示成一个函数，然后DEM格网的高程可表达为这个函数的积分平均表达式，利用已知GPS实测点与相邻格网建立平差模型，再利用最小二乘解算出格网的高程改正值和四个角点的高程值，接着把四个角点当做下一次迭代的已知高程值进行联合解算，从而实现提高 SRTM精度的目的。
附图说明
图1为DEM示意图。
图2为实施例1中平坦区的SRTM影像图。
图3为实施例1中半平坦区的SRTM影像图。
图4为实施例1中陡峭区的SRTM影像图。
图5为实施例1中已知GPS观测点对应的像素位置图。
图6为实施例1中平坦区SRTM高程与GPS高程直线拟合图。
图7为实施例1中平坦区算法改正高程与GPS高程直线拟合图。
图8为实施例1中半平坦区SRTM高程与GPS高程直线拟合图。
图9为实施例1中半平坦区算法改正高程与GPS高程直线拟合图。
图10实施例1中为陡峭区SRTM高程与GPS高程直线拟合图。
图11实施例1中为陡峭区算法改正高程与GPS高程直线拟合图。
具体实施方式
为了便于理解本发明，下文将结合较佳的实施例对本发明作更全面、细致地描述，但本发明的保护范围并不限于以下具体的实施例。
除非另有定义，下文中所使用的所有专业术语与本领域技术人员通常理解的含义相同。本文中所使用的专业术语只是为了描述具体实施例的目的，并不是旨在限制本发明的保护范围。
实施例1
为了证明该算法的可行性与有效性，选取中国某地区为研究区域，从网站http://earthexplorer.usgs.gov/上面下载该地区的SRTM数据，然后从该地区截取三个地理特征不同的区域作为实验数据，它们的大小均为64*64，三个区域的坡度分别在0-5度，5-20 度，20度以上，如图2所示，分别对应平坦区、半平坦区、陡峭区，地理特征极具代表性。该区域已知高程点的高程来自于CORS站记录的该区域大量的GPS流动站点的高程。首先对三个研究区域分别任意选择一个已知的GPS观测点，已知GPS观测点与研究区域的位置关系如图3所示，然后利用该观测点获得的GPS高程值采用点面融合积分平差算法改进对应研究区域的SRTM影像。对于研究区域已知的其他GPS高程值则作为精度评定样本，则三个研究区域对应的GPS校验点分别为48、71、37个。
为了客观地评估改进的DEM精度，考虑到系统偏差和标准差等精度估计会受到异常值和误差非正态分布的影响，因此这里考虑采用稳健统计方法，采用误差均值(MeanError，ME)、均方根误差(Root Mean Square Error，RMSE)、归一化绝对偏差中位数(Normalized Median， NMAD)和误差的标准差(Standard deviation of error，SDE)来进行精度评定。
Figure BDA0002098861780000061
Figure BDA0002098861780000062
NMAD＝1.4826×medianj(|Δhj-mΔh|)
Figure BDA0002098861780000063
其中，N是GPS精度评定样本的个数；Hi，Href分别是提出算法改进的GPS高程和GPS的高程；Δhi,mΔh分别是第i个精度评定点的高程差和高程差的中值；这些值都是精度的度量，越小表明改进算法效果越好。
为避免偶然性，对每一个研究区域将进行多次独立重复实验。对研究区域GPS精度评定样本点与改进高程值比对，得到精度评定如表1所示，表中的精度为多次实验的平均值。
表1三个研究区域精度评定
Figure BDA0002098861780000064
由表1可得，在平坦区、过渡区和陡峭区采用所提算法得到的DEM具有明显优势，与SRTM DEM 比较，改进的RMSE提高了1-2米，精度分别提高了19.5％、20.3％和11.6％，表明所提出算法能够有效改善SRTM DEM，特别是在平坦区域和和半平坦区域，但是在陡峭区域，因为受到地势的影响，融合的精度将受到一定程度的影响。为了更好地展示和分析所提出算法的优势，这里对精度评定GPS点的测量高程和采用所提出算法对这些点融合之后的高程之间的测量关系进行详细分析，分别列出GPS测量高程与相应SRTM的高程的直线拟合图，和GPS测量高程与提出算法改正高程的直线拟合图，如图6-图11所示。
图6、图8和图10为三个研究区域表明精度评定点的SRTM高程值与GPS观测高程之间的相关性。图7、图9和图11反映的是算法改正高程与GPS观测高程之间的相关性。这六个图显示的均是线性与正斜率的关系，表明两两变量之间朝着同一个方向发展。图7、图9和图11分别与图6、图8和图10比拥有更高的相关性，而且散点的分布更加具有一致性，说明提出的点面融合的平差算法改进DEM效果明显。
本实施例中，结果表明所提出点面融合的积分平差算法能够有效改善平坦区和和半平坦区的SRTM DEM，精度大约提高20％，但是在陡峭区域，因为受到地势的影响，融合的精度将受到一定程度的影响，精度大约能提高11％。

Claims (4)

1.一种GPS与SRTM点面融合的方法，其特征在于，包括以下步骤：
1)通过地面GPS观测，获得目标地区SRTM的DEM格网上任一点p的GPS高程值hp；
2)用测量平差方法把hp和p所在DEM格网的高程值
Figure FDA0003944809840000011
进行融合；得到p点所在格网的格网值和四个角点a、b、c、d新的高程值；作为第一层融合；
3)以a、b、c、d四个点为新的观测值，分别对其周围的格网进行融合，得到外围8格网融合后新的格网值及外围12个节点的新的高程；作为第二层融合；计算第二层8个格网平差融合后的新值与原格网值的差值，若各差值的均值<阀值，则计算停止；若各差值的均值>阀值的格网，则修正原DEM的高程值，直到各差值的均值<阀值；依次类推，得到全格网的高程值；当有格网的角点得到多个高程值时，取平均值。
2.根据权利要求1所述的GPS与SRTM点面融合的方法，其特征在于：所述DEM中每个格网的高程值
Figure FDA0003944809840000012
表示这个格网内所有点高程的平均值；
用h(x,y)表示任意一点(x,y)的高程，
Figure FDA0003944809840000013
表示格网(i,j)的格网值,那么
Figure FDA0003944809840000014
可以用h(x,y)在这个格网内的积分的平均值表示，即：
Figure FDA0003944809840000015
其中，S为格网面积。
3.根据权利要求2所述的GPS与SRTM点面融合的方法，其特征在于：所述DEM的各格网值是通过各种测绘手段获取的
Figure FDA0003944809840000016
的观测值hi,j；考虑观测误差，有：
Figure FDA0003944809840000017
p点的高程可表示为h(xp，yp)，考虑到观测误差的存在，p点高程的观测值可表示为：
hp+v＝h(xp，yp) (3) 。
4.根据权利要求1所述的GPS与SRTM点面融合的方法，其特征在于：将目标地区分成N个单元，每个单元独立的进行hp和p所在DEM格网的高程值
Figure FDA0003944809840000018
的融合。
CN201910528151.7A 2019-06-18 2019-06-18 一种gps与srtm点面融合的方法 Active CN110310370B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910528151.7A CN110310370B (zh) 2019-06-18 2019-06-18 一种gps与srtm点面融合的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910528151.7A CN110310370B (zh) 2019-06-18 2019-06-18 一种gps与srtm点面融合的方法

Publications (2)

Publication Number Publication Date
CN110310370A CN110310370A (zh) 2019-10-08
CN110310370B true CN110310370B (zh) 2023-01-10

Family

ID=68076617

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910528151.7A Active CN110310370B (zh) 2019-06-18 2019-06-18 一种gps与srtm点面融合的方法

Country Status (1)

Country Link
CN (1) CN110310370B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111595293A (zh) * 2020-05-29 2020-08-28 山东交通学院 一种基于多源监测数据融合的地表形变监测方法与系统
CN114398592B (zh) * 2021-12-27 2022-09-02 中国人民武装警察部队警官学院 一种基于异质格网高程线性分解模型的高程值解算方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006064051A1 (fr) * 2004-12-16 2006-06-22 Centre National D'etudes Spatiales Procede de traitement d'images mettant en oeuvre le georeferencement automatique d'images issues d'un couple d'images pris dans le meme plan focal
CN103823981A (zh) * 2014-02-28 2014-05-28 武汉大学 一种数字高程模型辅助的卫星影像区域网平差方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006064051A1 (fr) * 2004-12-16 2006-06-22 Centre National D'etudes Spatiales Procede de traitement d'images mettant en oeuvre le georeferencement automatique d'images issues d'un couple d'images pris dans le meme plan focal
CN103823981A (zh) * 2014-02-28 2014-05-28 武汉大学 一种数字高程模型辅助的卫星影像区域网平差方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
《Establishing Water Surface Area-Storage Capacity Relationship of Small Tanks Using SRTM and GPS》;V.Venkatesan等;《Energy Procedia》;20120316;第1167-1173页 *
《Improving the Accuracy of the SRTM Global DEM Using GPS data fusion and regression Model》;Mona Saad ElSayed等;《International Journal of Engineering Research》;20161110;第190-196页 *
SRTM约束的无地面控制立体影像区域网平差;周平等;《测绘学报》;20161115(第11期);第62-71页 *
多源DEM融合的高差拟合神经网络方法;沈焕锋等;《测绘学报》;20180615(第06期);第168-177页 *

Also Published As

Publication number Publication date
CN110310370A (zh) 2019-10-08

Similar Documents

Publication Publication Date Title
Elkhrachy Vertical accuracy assessment for SRTM and ASTER Digital Elevation Models: A case study of Najran city, Saudi Arabia
Nikolakopoulos et al. SRTM vs ASTER elevation products. Comparison for two regions in Crete, Greece
Catalão et al. Merging GPS and atmospherically corrected InSAR data to map 3-D terrain displacement velocity
CN103823981B (zh) 一种数字高程模型辅助的卫星影像区域网平差方法
CN102506824B (zh) 一种城市低空无人机系统生成数字正射影像图的方法
Pakoksung et al. Digital elevation models on accuracy validation and bias correction in vertical
CN109100719B (zh) 基于星载sar影像与光学影像的地形图联合测图方法
Tang et al. Verification of ZY-3 satellite imagery geometric accuracy without ground control points
CN111105496A (zh) 一种基于机载激光雷达点云数据的高精dem构建方法
Cao et al. Bundle adjustment of satellite images based on an equivalent geometric sensor model with digital elevation model
Chrysoulakis et al. Validation of ASTER GDEM for the Area of Greece
Agarwal et al. Accuracy assessment of the CARTOSAT DEM using robust statistical measures
CN110310370B (zh) 一种gps与srtm点面融合的方法
Gonçalves et al. Monitoring Local Shoreline Changes by Integrating UASs, Airborne LiDAR, Historical Images and Orthophotos.
De Oliveira et al. Assessment of radargrammetric DSMs from TerraSAR-X Stripmap images in a mountainous relief area of the Amazon region
Mukherjee et al. Accuracy of Cartosat-1 DEM and its derived attribute at multiple scale representation
Adavi et al. Analyzing different parameterization methods in GNSS tomography using the COST benchmark dataset
Isioye et al. Comparison and validation of ASTER-GDEM and SRTM elevation models over parts of Kaduna State, Nigeria
Zhou et al. DEM generation using Ziyuan-3 mapping satellite imagery without ground control points
Hassan et al. Cadastral mapping accuracy assessment using various surveying techniques and high-resolution satellites images
Santos Filho et al. Cartographic Accuracy Standard (CAS) of the digital terrain model of the digital and continuous cartographic base of the state of Amapá: case study in the city of Macapá
Sefercik et al. Point-based and model-based geolocation analysis of airborne laser scanning data
CN110322557B (zh) 一种结合gps与srtm融合的游动平差方法
Mo et al. Evaluation of AW3D30 Elevation Accuracy in China.
Rizeei et al. Extraction and accuracy assessment of DTMs derived from remotely sensed and field surveying approaches in GIS framework

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