CN110146924B - 基于水波初至偏振方位的海底地震仪位置及方位反演方法 - Google Patents

基于水波初至偏振方位的海底地震仪位置及方位反演方法 Download PDF

Info

Publication number
CN110146924B
CN110146924B CN201910598077.6A CN201910598077A CN110146924B CN 110146924 B CN110146924 B CN 110146924B CN 201910598077 A CN201910598077 A CN 201910598077A CN 110146924 B CN110146924 B CN 110146924B
Authority
CN
China
Prior art keywords
corrected
seismograph
azimuth
ocean bottom
orientation
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
CN201910598077.6A
Other languages
English (en)
Other versions
CN110146924A (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 University of Geosciences Beijing
Original Assignee
China University of Geosciences Beijing
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 University of Geosciences Beijing filed Critical China University of Geosciences Beijing
Priority to CN201910598077.6A priority Critical patent/CN110146924B/zh
Publication of CN110146924A publication Critical patent/CN110146924A/zh
Application granted granted Critical
Publication of CN110146924B publication Critical patent/CN110146924B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/38Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas
    • G01V1/3843Deployment of seismic devices, e.g. of streamers
    • G01V1/3852Deployment of seismic devices, e.g. of streamers to the seabed
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/10Aspects of acoustic signal generation or detection
    • G01V2210/14Signal detection
    • G01V2210/142Receiver location
    • G01V2210/1427Sea bed

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Oceanography (AREA)
  • Engineering & Computer Science (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种基于水波初至偏振方位的海底地震仪位置及方位反演方法,采用非线性反演方法,对于一个待修正的海底地震仪,利用水波初至偏振方位残差之差构造目标函数,并求取令目标函数最小的地震仪位置,在求得精确的地震仪位置之后,再以炮点方位与偏振方位之差求得地震仪的实际方位指向。尤其是以多对水波初至偏振方位残差之差的平方和作为目标函数,应用DE最优化方法求取使目标函数取最小值的地震仪位置的方案,在目标函数中消去了地震仪方位的影响,具有较快的收敛速度和较好的稳定性,相比线性方法具有较好的精度。

Description

基于水波初至偏振方位的海底地震仪位置及方位反演方法
技术领域
本发明涉及海底地震仪的位置及方法反演方法,具体涉及一种基于水波初至偏振方位的海底地震仪位置及方位反演方法。
背景技术
OBS海底地震仪在布放到海底之后,其方位指向通常是随机的,由于海流等因素的影响,其实际位置通常也会偏离预定位置,因此需要对其位置及方位进行反演处理。传统方法采用线性反演的方式,在一些特殊情况下,存在不收敛的情况。
发明内容
针对现有技术的不足,本发明旨在提供一种基于水波初至偏振方位的海底地震仪位置及方位反演方法,采用非线性反演方法,对于一个待修正的海底地震仪,以多对水波初至偏振方位残差之差的平方和作为目标函数,应用DE最优化方法求取使目标函数取最小值的地震仪位置,在目标函数中消去了地震仪方位的影响,具有较快的收敛速度和较好的稳定性。
为了实现上述技术目的,本发明采用如下技术方案:
基于水波初至偏振方位的海底地震仪位置及方位反演方法,具体为:
S1、根据待修正的海底地震仪的设计布放位置,利用该地震仪接收N个不同位置的炮点的三分量地震波形信号;
S2、对于一个待修正的海底地震仪,利用信号的水波初至偏振方位残差构造目标函数,并求取使目标函数取最小值的待修正的海底地震仪位置;
S3、利用步骤S1求得使目标函数取最小值的待修正的海底地震仪位置计算待修正的海底地震仪位置指向炮点方向的方位角,再以该待修正的海底地震仪位置指向炮点的方位角与炮点信号的水波初至偏振方位角之差求得待修正的海底地震仪的实际方位指向,并据此得到海底地震仪方位的修正量。
进一步地,步骤S2中,以多对信号的水波初至偏振方位残差之差的平方和作为目标函数,应用DE最优化方法求取使目标函数取最小值的待修正的海底地震仪位置。
进一步地,步骤S2的具体过程为:
设待修正的海底地震仪在海底的设计布放位置为(x,y),实际位置为(x+δx,y+δy),(δx,δy)为偏移量;
用多个炮点对的方位残差之差的平方和构造目标函数S2
Figure BDA0002117661210000021
Figure BDA0002117661210000022
Figure BDA0002117661210000023
别为待修正的海底地震仪位置指向第i号和第j号炮点的方位角的观测值,i≠j,表达式如下:
Figure BDA0002117661210000024
Figure BDA0002117661210000025
Figure BDA0002117661210000026
Figure BDA0002117661210000027
分别为待修正的海底地震仪位置指向第i号和第j号炮点方向的方位角的理论值,i≠j,表达式如下:
Figure BDA0002117661210000028
Figure BDA0002117661210000029
φ0表示待修正的海底地震仪的方位的修正量;
Figure BDA00021176612100000210
Figure BDA00021176612100000211
分别为根据第i号和第j号炮点的信号的三分量水波初至波形计算得到的偏振方位角,
Figure BDA00021176612100000212
Figure BDA00021176612100000213
分别表示待修正的海底地震仪在初始位置偏离(δx,δy)时指向第i号和第j号炮点的理论方位角;
由此可得:
Figure BDA00021176612100000214
采用DE非线性最优化方法求取使S2取最小值的偏移量(δx,δy)。
作为另一种实施方式,步骤S2中,设待修正的海底地震仪在海底的设计布放位置为(x,y),实际位置为(x+δx,y+δy),(δx,δy)为偏移量;
构造目标函数S1,求取使S1取最小值的(φ0,δx,δy):
Figure BDA0002117661210000031
其中,
Figure BDA0002117661210000032
为待修正的海底地震仪位置指向第i号炮点的方位角的观测值,表达式如下:
Figure BDA0002117661210000033
Figure BDA0002117661210000034
为待修正的海底地震仪位置指向第i号炮点方向的方位角的理论值,表达式如下:
Figure BDA0002117661210000035
S1是所有炮点的方位残差的平方和,
Figure BDA0002117661210000036
表示在初始位置偏离(δx,δy)时指向i号炮点的理论方位角,N为炮点的数量;φ0表示待修正的海底地震仪的方位的修正量;
Figure BDA0002117661210000037
为根据第i号炮点的信号的三分量水波初至波形计算得到的偏振方位角。
更进一步地,步骤S3具体为:
针对第i号炮点计算的方位修正量φ0为:
Figure BDA0002117661210000038
对于一个待修正的海底地震仪,针对每个炮点i,1≤i≤N都可计算得到一个方位修正量φ0,对所有结果求平均值,即得到最终的海底地震仪方位修正量φ'0
Figure BDA0002117661210000039
本发明的有益效果在于:采用非线性反演方法,对于一个待修正的海底地震仪,利用水波初至偏振方位残差之差构造目标函数,并求取令目标函数最小的地震仪位置,在求得精确的地震仪位置之后,再以炮点方位与偏振方位之差求得地震仪的实际方位指向,准确性较好。尤其是以多对水波初至偏振方位残差之差的平方和作为目标函数,应用DE最优化方法求取使目标函数取最小值的地震仪位置的方案,在目标函数中消去了地震仪方位的影响,具有较快的收敛速度和较好的稳定性。
附图说明
图1为本发明实施例3中实际数据采用第二种方法的迭代收敛情况示意图;
图2为本发明实施例3中实际数据采用第一种方法的迭代收敛情况示意图。
具体实施方式
以下将结合附图对本发明作进一步的描述,需要说明的是,本实施例以本技术方案为前提,给出了详细的实施方式和具体的操作过程,但本发明的保护范围并不限于本实施例。
实施例1
本实施例提供一种基于水波初至偏振方位的海底地震仪位置及方位反演方法,具体为:
S1、根据待修正的海底地震仪的设计布放位置,利用该地震仪接收N个不同位置的炮点的三分量地震波形信号
S2、设待修正的海底地震仪在海底的设计布放位置为(x,y),实际位置为(x+δx,y+δy),(δx,δy)为偏移量。
根据纵波的性质,纵波的偏振方向与传播方向相同,因此炮点放出的信号的水波初至偏振方位与待修正的海底地震仪位置指向炮点的方向一致,待修正的海底地震仪位置指向第i号炮点方向的观测值
Figure BDA0002117661210000041
可以表示为:
Figure BDA0002117661210000042
其中
Figure BDA0002117661210000043
为根据第i号炮点的三分量水波初至波形计算得到的偏振方位角,φ0表示待修正的海底地震仪的方位的修正量;待修正的海底地震仪位置指向第i号炮点方向的理论值
Figure BDA0002117661210000051
可由待修正的海底地震仪位置与炮点位置计算得到:
Figure BDA0002117661210000052
本实施例提供两种方式计算待修正的海底地震仪的方位的修正量φ0
1、作为一种实施方式,直接的方法是,构造目标函数S1,求取使S1取最小值的(φ0,δx,δy):
Figure BDA0002117661210000053
S1是所有炮点方位残差的平方和,
Figure BDA0002117661210000054
表示在初始位置偏离(δx,δy)时指向i号炮点的理论方位角,N为炮点的数量。
2、作为另一种实施方式,可以用多个炮点对的方位残差之差的平方和构造目标函数S2
Figure BDA0002117661210000055
将前述各式代入上式,得到:
Figure BDA0002117661210000056
φ0表示待修正的海底地震仪的方位的修正量;
Figure BDA0002117661210000057
Figure BDA0002117661210000058
分别为根据第i号和第j号炮点的信号的三分量水波初至波形计算得到的偏振方位角,
Figure BDA0002117661210000059
Figure BDA00021176612100000510
分别表示待修正的海底地震仪在初始位置偏离(δx,δy)时指向第i号和第j号炮点的理论方位角;
式(5)中,对于同一待修正的海底地震仪,方位修正量相同,因而被消掉,因此只含有两个未知数。
采用DE(deferential evolution)非线性最优化方法求取使S2取最小值的偏移量(δx,δy),针对第i号炮点计算的方位修正量φ0为:
Figure BDA00021176612100000511
对于一个待修正的海底地震仪,针对每个炮点i,1≤i≤N都可计算得到一个方位修正量φ0,对所有结果求平均值,即得到最终的海底地震仪方位修正量φ仩0
Figure BDA0002117661210000061
实施例2
本实施例旨在提供一种如实施例1所述方法所得反演结果的误差计算方法。
本实施例采用bootstrap方法给出反演结果的误差估计,具体做法是根据给定的实际初至偏振数据和方位角误差信息,按实际的炮点/台站位置生成一组理论观测数据,并将原有数据的残差随机加到理论数据中,然后对这些理论观测数据进行反演计算,以结果的统计方差统作为误差。
对于一个海底地震仪,假定有N个炮点数据参与计算,所涉及数据包括台站水平位置及方位(lat,lon,φ0),炮点坐标、方位角及残差(lati,loni,θi,ri),根据这组数据生成一组“新的理论数据”的具体步骤如下:
1.生成随机序列Ri,i=1,N。先对Ri序列赋初值Ri=i;
2.对Ri中的每个元素进行循环,每次循环生成[1,N]区间内的一个随机数j,然后Ri和Rj的值对调。循环完成,则随机序列Ri确定。实际程序中采用C语言中的系统调用函数srand()/rand()生成随机数,并以调用时的系统时间作为seed。
3.则新的一组观测数据为
Figure BDA0002117661210000062
重复上述步骤生成M组“观测数据”,对每组数据进行反演得到M个结果,这些结果的统计方差即是反演结果的误差。
实施例3
本实施例旨在对实施例1所述的方法进行进一步的检验。
为检验上述方法,本实施例分别进行了理论模型试验和实际数据的处理,所用数据来自2003年在我国南海某区域进行的OBS观测,这次观测有5台OBS台站。理论模型采用其中的obs4的实际观测系统,在理论计算得到的方位角中加上均值为0并且符合高斯分布的随机扰动,随机扰动的方差为6°,(φ0,δx,δy)各值的真实值均为0。表1是进行了1000次这种计算的统计结果,从表1中的结果可以看出,非线性反演方法的结果好于线性反演的结果,各反演量的平均值更接近真实值,方差也较小。
表2是5个台站实际数据的反演结果。从结果可以看出,两种方法得到的地震仪方位φ0非常接近,而地震仪的位置有些台站两种方法偏差较大,达到30米左右,但结合表1理论模型的结果,这个偏差也在正常范围内。
表1
Figure BDA0002117661210000071
表2实际数据的反演结果
Figure BDA0002117661210000072
Figure BDA0002117661210000081
理论模型实验和实际数据的结果表明,采用公式(3)(5)的两种方法得到的结果几乎完全相同,但迭代收敛情况公式(5)明显好于公式(3)的方法。如图1-2所示,图1-2是一个实际的obs台站数据分别采用公式(5)(3)的迭代收敛情况,公式(5)在迭代到600多次时即达到最优稳定状态,而采用公式(3)的方法,在迭代1000多次之后才达到稳定最优状态。
相关参考文献:
[1]Christeson,G.L.,K.D.McIntosh,T.H.Shipley,E.R.Flueh,and H.Goedde(1999),Structure of the costa rica convergent margin,offshore nicoyapeninsula,JOURNAL OF GEOPHYSICAL RESEARCH-SOLID EARTH,104(B11),25,443–25,468.
[2]JURKEVICS,A.(1988),Polarization analysis of 3-component arraydata,B.SEISMOL.SOC.AM.,78(5),1725–1743.
[3]Li,X.Y.,andJ.X.Yuan(1999),Geophone orientation and coupling inthree-component sea-floor data:a case study,GEOPHYS.PROSPECT.,47(6),995–1013.
[4]NAKAMURA,Y.,P.L.DONOHO,P.H.ROPER,and P.M.MCPHERSON(1987),Large-offset seismic surveying using ocean-bottom seismographs and air guns-instrumentation and field technique,GEOPHYSICS,52(12),1601–1611.
[5]PARK,J.,F.L.VERNON,andC.R.LINDBERG(1987),Frequency-dependentpolarization analysis of high-frequency seismograms,J GEOPHYS RES-SOLID,92(B12),12,664–12,674.
[6]Ruzek,B.,and M.Kvasnicka(2001),Differential evolution algorithm inthe earthquake hypocenter location,PURE APPL GEOPHYS,158(4),667–693.
[7]Ruzek,B.,P.Hrubcova,M.Novotny,A.Spicak,and O.Karousova(2007),Inversion of travel times obtained during active seismic refractionexperiments celebration 2000,alp 2002 and sudetes 2003,STUD GEOPHYS GEOD,51(1),141–164.
对于本领域的技术人员来说,可以根据以上的技术方案和构思,给出各种相应的改变和变形,而所有的这些改变和变形,都应该包括在本发明权利要求的保护范围之内。

Claims (2)

1.基于水波初至偏振方位的海底地震仪位置及方位反演方法,其特征在于,具体为:
S1、根据待修正的海底地震仪的设计布放位置,利用该地震仪接收N个不同位置的炮点的三分量地震波形信号;
S2、对于一个待修正的海底地震仪,利用信号的水波初至偏振方位残差构造目标函数,并求取使目标函数取最小值的待修正的海底地震仪位置;
S3、利用步骤S2求得使目标函数取最小值的待修正的海底地震仪位置计算待修正的海底地震仪位置指向炮点方向的方位角,再以该待修正的海底地震仪位置指向炮点的方位角与炮点信号的水波初至偏振方位角之差求得待修正的海底地震仪的实际方位指向,并据此得到海底地震仪方位的修正量;
其中,
步骤S2的具体过程为:
设待修正的海底地震仪在海底的设计布放位置为(x,y),实际位置为(x+δx,y+δy),(δx,δy)为偏移量;
用多个炮点对的方位残差之差的平方和构造目标函数S2
Figure FDA0002450917580000011
Figure FDA0002450917580000012
Figure FDA0002450917580000013
分别为待修正的海底地震仪位置指向第i号和第j号炮点的方位角的观测值,i≠j,表达式如下:
Figure FDA0002450917580000014
Figure FDA0002450917580000015
Figure FDA0002450917580000016
Figure FDA0002450917580000017
分别为待修正的海底地震仪位置指向第i号和第j号炮点方向的方位角的理论值,i≠j,表达式如下:
Figure FDA0002450917580000018
Figure FDA0002450917580000021
φ0表示待修正的海底地震仪的方位的修正量;
Figure FDA0002450917580000022
Figure FDA0002450917580000023
分别为根据第i号和第j号炮点的信号的三分量水波初至波形计算得到的偏振方位角,
Figure FDA0002450917580000024
Figure FDA0002450917580000025
分别表示待修正的海底地震仪在初始位置偏离(δx,δy)时指向第i号和第j号炮点的理论方位角;
由此可得:
Figure FDA0002450917580000026
采用DE非线性最优化方法求取使S2取最小值的偏移量(δx,δy)。
2.根据权利要求1所述的基于水波初至偏振方位的海底地震仪位置及方位反演方法,其特征在于,
步骤S3具体为:
针对第i号炮点计算的方位修正量φ0为:
Figure FDA0002450917580000027
对于一个待修正的海底地震仪,针对每个炮点i,1≤i≤N都可计算得到一个方位修正量φ0,对所有结果求平均值,即得到最终的海底地震仪方位修正量φ'0
Figure FDA0002450917580000028
CN201910598077.6A 2019-07-03 2019-07-03 基于水波初至偏振方位的海底地震仪位置及方位反演方法 Active CN110146924B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910598077.6A CN110146924B (zh) 2019-07-03 2019-07-03 基于水波初至偏振方位的海底地震仪位置及方位反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910598077.6A CN110146924B (zh) 2019-07-03 2019-07-03 基于水波初至偏振方位的海底地震仪位置及方位反演方法

Publications (2)

Publication Number Publication Date
CN110146924A CN110146924A (zh) 2019-08-20
CN110146924B true CN110146924B (zh) 2020-05-26

Family

ID=67596916

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910598077.6A Active CN110146924B (zh) 2019-07-03 2019-07-03 基于水波初至偏振方位的海底地震仪位置及方位反演方法

Country Status (1)

Country Link
CN (1) CN110146924B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110824551B (zh) * 2019-11-14 2020-12-11 自然资源部第二海洋研究所 一种基于声学测距和多波束地形的obs精确定位方法
CN111257941B (zh) * 2020-02-18 2020-09-04 中国科学院地质与地球物理研究所 一种组合式海底地震仪方位角自动识别装置及方法
CN111352160B (zh) * 2020-03-19 2020-11-10 中国科学院地质与地球物理研究所 一种海底地震仪自动重定位装置及方法
CN115308800B (zh) * 2022-08-29 2024-07-12 南方海洋科学与工程广东省实验室(广州) 一种利用海底反射波走时和地形数据定位海底地震仪的方法及处理终端

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106154334A (zh) * 2015-04-13 2016-11-23 中石化石油工程地球物理有限公司胜利分公司 基于网格搜索的井下微地震事件实时反演定位方法
CN106249295A (zh) * 2015-06-15 2016-12-21 中国石油化工股份有限公司 一种井中微地震p、s波联合快速定位方法及系统
WO2017062322A1 (en) * 2015-10-07 2017-04-13 Schlumberger Technology Corporation Seismic sensor orientation

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106154334A (zh) * 2015-04-13 2016-11-23 中石化石油工程地球物理有限公司胜利分公司 基于网格搜索的井下微地震事件实时反演定位方法
CN106249295A (zh) * 2015-06-15 2016-12-21 中国石油化工股份有限公司 一种井中微地震p、s波联合快速定位方法及系统
WO2017062322A1 (en) * 2015-10-07 2017-04-13 Schlumberger Technology Corporation Seismic sensor orientation

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Determination of New Zealand Ocean Bottom Seismometer Orientation via Rayleigh-Wave Polarization;J. C. Stachnik et al.;《Seismological Research Letters》;20120831;第83卷(第4期);第704-712页 *
Ocean-Bottom Seismometer Instrument Orientations via Automated Rayleigh-Wave Arrival-Angle Measurements;Adrian K. Doran et al.;《Bulletin of the Seismological Society of America》;20170430;第107卷(第2期);第1-17页 *

Also Published As

Publication number Publication date
CN110146924A (zh) 2019-08-20

Similar Documents

Publication Publication Date Title
CN110146924B (zh) 基于水波初至偏振方位的海底地震仪位置及方位反演方法
CN108614268B (zh) 低空高速飞行目标的声学跟踪方法
CN104198992A (zh) 基于多径时延结构压缩感知的水声目标被动定位方法
CN109444956B (zh) 三维起伏观测面地震斜率层析成像方法
CN109765616B (zh) 一种保幅波场延拓校正方法及系统
CN108984818A (zh) 固定翼时间域航空电磁数据拟三维空间约束整体反演方法
CN110275139B (zh) 一种基于旋转式基元复用的超短基线定位系统及方法
CN108983169B (zh) 一种基于数字高程模型的米波雷达地形修正方法
CN103925904A (zh) 一种基于对称测线的超短基线安装角度偏差无偏估计方法
CN112099082B (zh) 一种共面元共方位角道集的地震回折波走时反演方法
CN114563778B (zh) 基于时延冗余测量的阵形校准方法
CN113740915B (zh) 一种球坐标系重力和接收函数联合反演地壳结构参数的方法
CN113960532B (zh) 一种基于假想源的二次定位计算的微地震定位方法
CN106054252B (zh) 一种叠前时间偏移的方法及装置
Lü et al. Application of GPS technology to build a mine-subsidence observation station
CN109212598B (zh) 基于直达波反演的三维空间二次定位方法
CN109633527A (zh) 基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法
CN112378376B (zh) 一种基于传感阵列和测斜仪的海底变形联合监测方法
Kuang et al. Robust constrained Kalman filter algorithm considering time registration for GNSS/acoustic joint positioning
CN110749861B (zh) 一种基于多假设深度的水底固定目标的三维定位方法
Liu et al. A new method for OBS relocation using direct water-wave arrival times from a shooting line and accurate bathymetric data
CN112147695B (zh) 一种海底节点检波器水下姿态定向方法
CN100349003C (zh) 一种消除起伏地表成像误差的方法
CN117368849A (zh) 一种基于炮点差分的全息二次定位方法及其应用
Sun et al. A New Angle-Calibration Method for Precise Ultra-Short Baseline Underwater Positioning.

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