CN112130151A - 一种圆弧合成孔径地基雷达坐标投影快速计算方法 - Google Patents

一种圆弧合成孔径地基雷达坐标投影快速计算方法 Download PDF

Info

Publication number
CN112130151A
CN112130151A CN202011110192.3A CN202011110192A CN112130151A CN 112130151 A CN112130151 A CN 112130151A CN 202011110192 A CN202011110192 A CN 202011110192A CN 112130151 A CN112130151 A CN 112130151A
Authority
CN
China
Prior art keywords
radar
projection
polar coordinate
monitoring
coordinate system
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
Application number
CN202011110192.3A
Other languages
English (en)
Other versions
CN112130151B (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.)
Chinese Nonferrous Metal Survey And Design Institute Of Changsha Co ltd
Original Assignee
Chinese Nonferrous Metal Survey And Design Institute Of Changsha Co ltd
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 Chinese Nonferrous Metal Survey And Design Institute Of Changsha Co ltd filed Critical Chinese Nonferrous Metal Survey And Design Institute Of Changsha Co ltd
Priority to CN202011110192.3A priority Critical patent/CN112130151B/zh
Publication of CN112130151A publication Critical patent/CN112130151A/zh
Application granted granted Critical
Publication of CN112130151B publication Critical patent/CN112130151B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明提供了一种圆弧合成孔径地基雷达坐标投影快速计算方法。包括以下步骤:获取雷达照射场景在以雷达为中心的极坐标系下的坐标值(r,θ)序列以及地形三维点云数据;计算地形点云在此坐标系下的投影极坐标(rd,θd);选取一定数量的地形点云的投影极坐标值,对于每个地形点云的投影极坐标值,若||rd‑ri||<Δr,||θd‑θi||<Δθ,则认为该点云为照射场景在雷达极坐标系下的某一个坐标值为(ri,θi)的测量网格所覆盖的地形点云。本发明能够将照射场景在雷达极坐标系下坐标投影到地形上,能够在三维系统中直观的展示雷达监测范围及形变值。同时能够将每个监测网格中覆盖的地形点云筛选出来,在进行三维模型切片时能够在不同缩放级别中都展示雷达监测数据。

Description

一种圆弧合成孔径地基雷达坐标投影快速计算方法
技术领域
本发明专利涉及测绘及软件领域,特别地,涉及一种圆弧合成孔径地基雷达坐标投影快速计算方法。
背景技术
边坡地质灾害预警一直是矿山监测、地质勘测、应急救援等应用的关键任务之一,而相对于位移计、北斗差分监测设备、水准仪等常规监测技术手段,采用地基雷达进行远距离遥感监测的方式不受雨、雪、雾等气候影响和强光、夜晚等光照影响,测量精度达亚毫米级,可对大范围场景进行全面和快速测量,综合使用成本低,是一种极有前景的边坡形变监测技术手段。
圆弧合成孔径雷达在距离方向形成圆形合成孔径,可以对目标区域360°全方位监测。其成像结果为扇形区域网格,距离向和角度向网格个数依赖于雷达角度分辨率、距离分辨率以及监测范围大小。以极坐标形式表示,其中角度为目标到雷达合成孔径中心所形成的线在达二维坐标平面下偏离雷达起始位置的角度,距离为目标到雷达合成孔径中心的距离。这种极坐标表示的图像不直观,很难判别是哪个区域内发生形变,因此需要将极坐标投影到监测区域的地形上。
发明专利内容
本发明目的在于提供一种圆弧合成孔径地基雷达坐标投影快速计算方法,能够实现雷达监测结果与地形叠加,大幅度降低计算量,提高计算速度。
为实现上述目的,本发明提供了一种圆弧合成孔径地基雷达坐标投影快速计算方法,包括以下步骤:
S1、获取包含监测区域在内的整个项目范围内地形三维点云模型,其坐标集合为P(Xd,Yd,Hd);
S2、获取雷达一次扫描数据,计算得到照射区域内所有监测网格在以雷达转轴中心为极坐标系中心的极坐标值集合R(rij,θij)(1≤i≤M,1≤j≤N),极坐标系方位角起算位置为雷达启动时初始方位角;
S3、在UTM坐标系下,计算点云中每个点在所述雷达成像平面上的投影极坐标,计算方法如下:
Figure BDA0002728338560000021
Figure BDA0002728338560000022
其中,方位角θ方向定义为Y轴正方向为起始方位,顺时针为正;
将计算的每个点的Xd,Yd,Hd,rd,θd存入数据库表T1中;
S4、从集合R中取出一组方位角相同的极坐标值集合K(rij,θij),集合K中下标j相等,其在UTM坐标系下方位角θ'=θij0
S5、在表T1中,筛选出满足||θi-θ'||<Δθ的集合M(Xd,Yd,Hd,rd,θd);
S6、对于集合K中的每一个坐标(rij,θ'),若||θd-θ'||<Δθ,||rd–rij||<Δr,Δθ,Δr为设置的阈值,将数组(i,j,Xd,Yd,Hd)作为一条数据存在数据库表T2中;
S7、重复S4-S7步骤,直到计算完集合R中所有数据;
S8、数据库表T2中,具有相同的i、j值的(Xd,Yd,Hd)集合即为第i行第j列监测单元格所投影覆盖的地形点云。
进一步地,采用圆弧合成孔径雷达,雷达照射场景所覆盖的地形区域为监测范围,呈扇形格网,每一个监测单元为扇形方格,其位置表示为以雷达转轴中心为极坐标系中心的极坐标值(r,θ),其中θ为相对于雷达初始位置偏转的角度,r为距离;r,θ是以弧形方格内雷达回波信号的平均值计算;整个监测范围有M×N个监测单元,M为距离方向个数,N为角度方向个数。
进一步地,包括照射场景在内的整个项目区域,通过无人机航拍或者三维激光扫描仪获取地形三维模型,将模型导出为点云,其坐标为UTM投影下坐标。
进一步地,θ0为初始时雷达摆臂在UTM平面坐标系下方位角,采用全站仪或者RTK方式测量得到;X0、Y0、H0为雷达转轴中心在UTM平面坐标系下坐标。
进一步地,数据库表T1和表T2均采用sqlite数据库。
进一步地,Δθ,Δr为设置的阈值,阈值大小分别选择雷达角度分辨的一半和距离分辨的一半。
应用本发明的技术方案,具有以下有益效果:
(1)本发明能够将照射场景在雷达极坐标系下坐标投影到地形上,实现在三维系统中展示时,能够直观的展示雷达监测范围及形变值。
(2)采用本发明方法,能够快实现雷达监测区域到地形的投影,减少了计算量,计算速度大大提高,而且能够将每个监测网格中覆盖的地形点云筛选出来,在进行三维模型切片时能够在不同缩放级别中都展示雷达监测数据。
附图说明
构成本申请的一部分的附图用来提供对本发明专利的进一步理解,本发明专利的示意性实施例及其说明用于解释本发明专利,并不构成对本发明专利的不当限定。在附图中:
图1是本发明专利雷达照射区域示意图;
其中,1、照射场景,2、监测单元。
具体实施方式
以下结合附图对本发明专利的实施例进行详细说明,但是本发明专利可以根据权利要求限定和覆盖的多种不同方式实施。
实施例1:
为了实现雷达监测结果与地形叠加,本发明提出了一种圆弧合成孔径地基雷达坐标投影快速计算方法。如图1所示,本计算方法,采用圆弧合成孔径雷达,雷达照射场景1所覆盖的地形区域为监测范围,呈扇形格网,每一个监测单元2为扇形方格,其位置表示为以雷达转轴中心为极坐标系中心的极坐标值(r,θ),其中θ为相对于雷达初始位置偏转的角度,r为距离。r,θ是以弧形方格内雷达回波信号的平均值计算。整个监测范围有M×N个监测单元,M为距离方向个数,N为角度方向个数。
包括照射场景在内的整个项目区域,通过无人机航拍或者三维激光扫描仪获取地形三维模型,将模型导出为点云,其坐标为UTM投影下坐标,至少包括X、Y、H三个方向的坐标。
需要计算初始时雷达摆臂在UTM平面坐标系下方位角θ0,采用全站仪或者RTK(载波相位差分技术)方式测量;需要计算雷达转轴中心在UTM平面坐标系下坐标X0、Y0、H0
本发明方法包括以下步骤:
S1、获取包含监测区域在内的整个项目范围内地形三维点云模型,其坐标集合为P(Xd,Yd,Hd);
S2、获取雷达一次扫描数据,计算得到照射区域内所有监测网格在以雷达转轴中心为极坐标系中心的极坐标值集合R(rij,θij)(1≤i≤M,1≤j≤N),极坐标系方位角起算位置为雷达启动时初始方位角。
S3、在UTM坐标系下,计算点云中每个点在所述雷达成像平面上的投影极坐标,计算方法如下:
Figure BDA0002728338560000041
Figure BDA0002728338560000042
其中,方位角θ方向定义为Y轴正方向为起始方位,顺时针为正。
将计算的每个点的Xd,Yd,Hd,rd,θd存入数据库表T1中。
S4、从集合R中取出一组方位角相同的极坐标值集合K(rij,θij)(其中j相等),其在UTM坐标系下方位角θ'=θij0
S5、在表T1中,筛选出满足||θi-θ'||<Δθ的集合M(Xd,Yd,Hd,rd,θd)。
S6、对于集合K中的每一个坐标(rij,θ'),若||θd-θ'||<Δθ,||rd–rij||<Δr,Δθ,Δr为设置的阈值,将数组(i,j,Xd,Yd,Hd)作为一条数据存在数据库表T2中。
S7、重复S4-S7步骤,直到计算完集合R中所有数据。
S8、数据库表T2中,具有相同的i、j值的(Xd,Yd,Hd)集合即为第i行第j列监测单元格所投影覆盖的地形点云。
所述数据库表T1和表T2均采用sqlite数据库。
所述Δθ,Δr为设置的阈值,一般选择雷达角度分辨和距离分辨的一半。
选取实例,雷达坐标为(498950.45,5573596.81,786.75),初始方位角为206.9615°。雷达监测单元数,角度方向360个,距离方向1187个,雷达角度分辨率0.00418rad,距离分辨率0.5m,Δθ设置为0.00209rad,Δr设置为0.5m。整个项目范围内三维点云模型包含点云个数约为143万。
采用本发明所述方法,计算时间约为10分钟;而采用常规逐点计算,计算时间>12小时。采用本发明方法,能够快实现雷达监测区域到地形的投影,减少了计算量,计算速度大大提高。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (6)

1.一种圆弧合成孔径地基雷达坐标投影快速计算方法,其特征在于,包括以下步骤:
S1、获取包含监测区域在内的整个项目范围内地形三维点云模型,其坐标集合为P(Xd,Yd,Hd);
S2、获取雷达一次扫描数据,计算得到照射区域内所有监测网格在以雷达转轴中心为极坐标系中心的极坐标值集合R(rij,θij)(1≤i≤M,1≤j≤N),极坐标系方位角起算位置为雷达启动时初始方位角;
S3、在UTM坐标系下,计算点云中每个点在所述雷达成像平面上的投影极坐标,计算方法如下:
Figure FDA0002728338550000012
其中,方位角θ方向定义为Y轴正方向为起始方位,顺时针为正;
将计算的每个点的Xd,Yd,Hd,rd,θd存入数据库表T1中;
S4、从集合R中取出一组方位角相同的极坐标值集合K(rij,θij),集合K中下标j相等,其在UTM坐标系下方位角θ'=θij0
S5、在表T1中,筛选出满足||θi-θ'||<Δθ的集合M(Xd,Yd,Hd,rd,θd);
S6、对于集合K中的每一个坐标(rij,θ'),若||θd-θ'||<Δθ,||rd–rij||<Δr,Δθ,Δr为设置的阈值,将数组(i,j,Xd,Yd,Hd)作为一条数据存在数据库表T2中;
S7、重复S4-S7步骤,直到计算完集合R中所有数据;
S8、数据库表T2中,具有相同的i、j值的(Xd,Yd,Hd)集合即为第i行第j列监测单元格所投影覆盖的地形点云。
2.根据权利要求1所述的一种圆弧合成孔径地基雷达坐标投影快速计算方法,其特征在于:采用圆弧合成孔径雷达,雷达照射场景所覆盖的地形区域为监测范围,呈扇形格网,每一个监测单元为扇形方格,其位置表示为以雷达转轴中心为极坐标系中心的极坐标值(r,θ),其中θ为相对于雷达初始位置偏转的角度,r为距离;r,θ是以弧形方格内雷达回波信号的平均值计算;整个监测范围有M×N个监测单元,M为距离方向个数,N为角度方向个数。
3.根据权利要求1所述的一种圆弧合成孔径地基雷达坐标投影快速计算方法,其特征在于:包括照射场景在内的整个项目区域,通过无人机航拍或者三维激光扫描仪获取地形三维模型,将模型导出为点云,其坐标为UTM投影下坐标。
4.根据权利要求1所述的一种圆弧合成孔径地基雷达坐标投影快速计算方法,其特征在于:θ0为初始时雷达摆臂在UTM平面坐标系下方位角,采用全站仪或者RTK方式测量得到;X0、Y0、H0为雷达转轴中心在UTM平面坐标系下坐标。
5.根据权利要求4述的一种圆弧合成孔径地基雷达坐标投影快速计算方法,其特征在于:数据库表T1和表T2均采用sqlite数据库。
6.根据权利要求5述的一种圆弧合成孔径地基雷达坐标投影快速计算方法,其特征在于:Δθ,Δr为设置的阈值,阈值大小一般分别选择雷达角度分辨的一半和距离分辨的一半。
CN202011110192.3A 2020-10-16 2020-10-16 一种圆弧合成孔径地基雷达坐标投影快速计算方法 Active CN112130151B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011110192.3A CN112130151B (zh) 2020-10-16 2020-10-16 一种圆弧合成孔径地基雷达坐标投影快速计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011110192.3A CN112130151B (zh) 2020-10-16 2020-10-16 一种圆弧合成孔径地基雷达坐标投影快速计算方法

Publications (2)

Publication Number Publication Date
CN112130151A true CN112130151A (zh) 2020-12-25
CN112130151B CN112130151B (zh) 2022-07-08

Family

ID=73853211

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011110192.3A Active CN112130151B (zh) 2020-10-16 2020-10-16 一种圆弧合成孔径地基雷达坐标投影快速计算方法

Country Status (1)

Country Link
CN (1) CN112130151B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113051304A (zh) * 2021-04-02 2021-06-29 中国有色金属长沙勘察设计研究院有限公司 一种雷达监测数据与三维点云融合的计算方法
CN113932703A (zh) * 2021-11-09 2022-01-14 中国有色金属长沙勘察设计研究院有限公司 一种形变监测雷达区域数据处理方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102096713A (zh) * 2011-01-29 2011-06-15 广州都市圈网络科技有限公司 一种基于网格化二三维地图匹配方法及系统
CN102645209A (zh) * 2012-04-24 2012-08-22 长江勘测规划设计研究有限责任公司 机载LiDAR点云和高分辨率影像进行空间点的联合定位方法
CN105931234A (zh) * 2016-04-19 2016-09-07 东北林业大学 一种地面三维激光扫描点云与影像融合及配准的方法
CN106530380A (zh) * 2016-09-20 2017-03-22 长安大学 一种基于三维激光雷达的地面点云分割方法
CN107316325A (zh) * 2017-06-07 2017-11-03 华南理工大学 一种基于图像配准的机载激光点云与影像配准融合方法
CN109345574A (zh) * 2018-08-31 2019-02-15 西安电子科技大学 基于语义点云配准的激光雷达三维建图方法
CN110244302A (zh) * 2019-07-05 2019-09-17 苏州科技大学 地基合成孔径雷达影像像元坐标三维变换方法
CN110764111A (zh) * 2019-11-15 2020-02-07 深圳市镭神智能系统有限公司 雷达坐标与大地坐标的转换方法、装置、系统及介质

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102096713A (zh) * 2011-01-29 2011-06-15 广州都市圈网络科技有限公司 一种基于网格化二三维地图匹配方法及系统
CN102645209A (zh) * 2012-04-24 2012-08-22 长江勘测规划设计研究有限责任公司 机载LiDAR点云和高分辨率影像进行空间点的联合定位方法
CN105931234A (zh) * 2016-04-19 2016-09-07 东北林业大学 一种地面三维激光扫描点云与影像融合及配准的方法
CN106530380A (zh) * 2016-09-20 2017-03-22 长安大学 一种基于三维激光雷达的地面点云分割方法
CN107316325A (zh) * 2017-06-07 2017-11-03 华南理工大学 一种基于图像配准的机载激光点云与影像配准融合方法
CN109345574A (zh) * 2018-08-31 2019-02-15 西安电子科技大学 基于语义点云配准的激光雷达三维建图方法
CN110244302A (zh) * 2019-07-05 2019-09-17 苏州科技大学 地基合成孔径雷达影像像元坐标三维变换方法
CN110764111A (zh) * 2019-11-15 2020-02-07 深圳市镭神智能系统有限公司 雷达坐标与大地坐标的转换方法、装置、系统及介质

Non-Patent Citations (9)

* Cited by examiner, † Cited by third party
Title
杨俊等: "地基SAR图像与地形数据的几何映射三维匹配方法", 《中国科学院大学学报》, no. 03, 15 May 2015 (2015-05-15) *
林等: "圆迹SAR极坐标格式算法研究", 《电子与信息学报》 *
林等: "圆迹SAR极坐标格式算法研究", 《电子与信息学报》, no. 12, 15 December 2010 (2010-12-15) *
王鹏等: "GB-SAR影像坐标到三维地形坐标转换方法", 《长江科学院院报》 *
王鹏等: "GB-SAR影像坐标到三维地形坐标转换方法", 《长江科学院院报》, no. 06, 15 June 2018 (2018-06-15) *
赵杰等: "网格坐标转换为雷达测量坐标的工程化算法", 《火力与指挥控制》 *
赵杰等: "网格坐标转换为雷达测量坐标的工程化算法", 《火力与指挥控制》, no. 07, 15 July 2007 (2007-07-15) *
陈乐平等: "圆周合成孔径雷达的快速时域成像算法", 《国防科技大学学报》, no. 02, 28 April 2018 (2018-04-28) *
高德章: "大地坐标系与投影坐标系", 《物探化探计算技术》, no. 01, 15 January 2011 (2011-01-15) *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113051304A (zh) * 2021-04-02 2021-06-29 中国有色金属长沙勘察设计研究院有限公司 一种雷达监测数据与三维点云融合的计算方法
CN113051304B (zh) * 2021-04-02 2022-06-24 中国有色金属长沙勘察设计研究院有限公司 一种雷达监测数据与三维点云融合的计算方法
CN113932703A (zh) * 2021-11-09 2022-01-14 中国有色金属长沙勘察设计研究院有限公司 一种形变监测雷达区域数据处理方法

Also Published As

Publication number Publication date
CN112130151B (zh) 2022-07-08

Similar Documents

Publication Publication Date Title
GREJNER‐BRZEZINSKA Direct exterior orientation of airborne imagery with GPS/INS system: Performance analysis
Li Mobile mapping: An emerging technology for spatial data acquisition
US9927513B2 (en) Method for determining the geographic coordinates of pixels in SAR images
CN108983248A (zh) 一种基于3d激光雷达及v2x的网联车定位方法
Nagihara et al. Use of a three‐dimensional laser scanner to digitally capture the topography of sand dunes in high spatial resolution
US10228456B2 (en) Determination of the position of a vehicle on or above a planet surface
CN112130151B (zh) 一种圆弧合成孔径地基雷达坐标投影快速计算方法
CN112184786B (zh) 一种基于合成视觉的目标定位方法
CN111522005A (zh) 一种形变监测与地形重构方法
CN103644918A (zh) 卫星对月探测数据定位处理方法
CN109631863A (zh) 一种空地结合的潮间带一体化测绘方法
CN111398980A (zh) 一种机载LiDAR数据处理的方法及装置
CN114721436A (zh) 一种面向无人机载高光谱成像系统的航线自动规划方法
CN114689015A (zh) 一种提高光学卫星立体影像dsm高程精度的方法
US8577608B2 (en) Method of obtaining a local terrain elevation database from detector means on board a vehicle
CN110160503B (zh) 一种顾及高程的无人机景观匹配定位方法
Yu et al. Automatic extrinsic self-calibration of mobile LiDAR systems based on planar and spherical features
Liu et al. Correction of positional errors and geometric distortions in topographic maps and DEMs using a rigorous SAR simulation technique
Li et al. A study of the potential attainable geometric accuracy of IKONOS satellite imagery
Li et al. Photogrammetric processing of high-resolution airborne and satellite linear array stereo images for mapping applications
AU2006300903A1 (en) Terrain mapping
Yousif et al. Accuracy enhancement of terrestrial mobile lidar data using theory of assimilation
CN113483739B (zh) 海上目标位置测量方法
US11428802B2 (en) Localization using particle filtering and image registration of radar against elevation datasets
Li Coastline mapping and change detection using one-meter resolution satellite imagery

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