CN113917526B - 基于非分裂完全匹配层吸收边界的正演模拟方法 - Google Patents

基于非分裂完全匹配层吸收边界的正演模拟方法 Download PDF

Info

Publication number
CN113917526B
CN113917526B CN202010661155.5A CN202010661155A CN113917526B CN 113917526 B CN113917526 B CN 113917526B CN 202010661155 A CN202010661155 A CN 202010661155A CN 113917526 B CN113917526 B CN 113917526B
Authority
CN
China
Prior art keywords
matching layer
split
iteration
equation
absorption boundary
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
CN202010661155.5A
Other languages
English (en)
Other versions
CN113917526A (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 Petroleum and Chemical Corp
Geophysical Research Institute of Sinopec Shengli Oilfield Co
Original Assignee
China Petroleum and Chemical Corp
Geophysical Research Institute of Sinopec Shengli Oilfield Co
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 Petroleum and Chemical Corp, Geophysical Research Institute of Sinopec Shengli Oilfield Co filed Critical China Petroleum and Chemical Corp
Priority to CN202010661155.5A priority Critical patent/CN113917526B/zh
Publication of CN113917526A publication Critical patent/CN113917526A/zh
Application granted granted Critical
Publication of CN113917526B publication Critical patent/CN113917526B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282Application of seismic models, synthetic seismograms
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/66Subsurface modeling
    • 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
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

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

Abstract

本发明提供一种基于非分裂完全匹配层吸收边界的正演模拟方法,该基于非分裂完全匹配层吸收边界的正演模拟方法包括:步骤1,建立二阶弹性波动方程;步骤2,利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界;步骤3,实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟。该基于非分裂完全匹配层吸收边界的正演模拟方法与常规的分裂完全匹配层吸收边界相比,非分裂完全匹配层吸收边界的吸收效果更好,能够吸收高角度入射波和掠射波,计算速度更快、所占内存更小,降低了算法对计算机内存的需求,更易满足实际生产需求。

Description

基于非分裂完全匹配层吸收边界的正演模拟方法
技术领域
本发明涉及地球物理勘查技术领域,特别是涉及到一种基于非分裂卷积完全匹配层弹性波数值模拟方法。
背景技术
波动方程法是应对日益复杂的地质构造地震波传播数值模拟的最重要方法。由于无法构建无限区域的模型,使用波动方程法进行地震波正演数值模拟时,一般需要用有限的区域模拟无限区域波传播过程,即添加人工边界。为了消除地震波传播到人工边界引起的虚假反射,就需要设置吸收边界条件,使得地震波可以无反射的穿过人工边界。因此,发展形成有效而稳定的吸收边界条件一直是地震波数值模拟的一项重要研究内容。完全匹配层吸收边界是在模型区域外人为增加了一层吸收介质,当波传播到该介质中时,波的能量将被介质以指数衰减的速度迅速吸收,直至消失,而不会反弹回模型区域。理论上,通过合理的设置该吸收层的吸收系数,无论什么频率的波以任意角度从模型区域进入吸收介质层时,都不产生反射。与传统的吸收边界条件相比,完全匹配层具有更优秀的吸收效果。目前应用较为广泛的完全匹配层吸收边界为分裂完全匹配层(SPML)吸收边界。但在大角度入射的情况下,SPML吸收边界的吸收效果并不理想,虚假反射将会更加明显,对薄层介质、震源靠近模型离散的网格边缘或者大偏移距模型的数值模拟效果不理想,此外,SPML边界无法吸收掠射波,并且正演时对内存需求较大。在分裂完全匹配层(SPML)吸收边界的基础上发展起来的卷积完全匹配层(convolution perfectly matched layer,CPML)吸收边界能够吸收大角度入射波和掠射波,但目前构建的卷积完全匹配层吸收边界多适用于一阶波动方程,难以直接应用于二阶波动方程。
为解决以上技术问题,本发明提出一种基于非分裂完全匹配层吸收边界正演模拟方法。
发明内容
本发明的目的是提供一种适用于二阶弹性波方程的非分裂卷积完全匹配层吸收边界条件的数值模拟方法。
本发明的目的可通过如下技术措施来实现:基于非分裂完全匹配层吸收边界的正演模拟方法,该基于非分裂完全匹配层吸收边界的正演模拟方法包括:步骤1,建立二阶弹性波动方程;步骤2,利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界;步骤3,实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟。
本发明的目的还可通过如下技术措施来实现:
在步骤1中,通过柯西方程、几何方程以及纳维尔方程推导出二阶弹性波波动方程:
其中λ和μ为介质的拉梅常数;ρ表示介质密度;uz为波场垂直分量,ux为波场水平分量,t表示地震波传播时间,x,z分别表示水平方向和垂直方向;
在步骤2中,采用复坐标变换,忽略衰减因子空变影响计算出二阶弹性波方程的空间二阶偏导算子:
其中:δ(t)和H(t)分别代表狄拉克函数以及单位阶跃函数,t为时间,dx、dz为衰减透射波的衰减参数;αx、αz为滤波的频率参数;为x方向在复频域的拉伸变量,/>为x与z方向在复频域的拉伸变量,kz,kz为衰减渐消波的尺度因子。
在步骤2中,利用迭代法代替卷积运算的思想,将二阶偏导算子中的时域卷积运算转换成迭代运算计算出时间域基于非分裂完全匹配层吸收边界的二阶弹性波方程。
在步骤3中,由连续方程的时间二阶、空间十阶有限差分计算基于非分裂完全匹配层吸收边界二阶弹性波方程的离散方程。
在步骤3中,对均匀介质模型和复杂层状模型分别开展分裂完全匹配层和非分裂完全匹配层吸收边界二阶精度弹性波方程正演模拟,通过成图验证非分裂完全匹配层吸收边界具有更好的吸收效果、更快的计算速度和对计算机内存较低的需求。
本发明的基于非分裂完全匹配层吸收边界的正演模拟方法,由连续方程的时间二阶、空间十阶有限差分计算基于非分裂完全匹配层吸收边界二阶弹性波方程的离散方程,基于卷积完全匹配层(CPML)吸收边界利用迭代运算取代卷积运算的思想推导出非分裂完全匹配层吸收边界(NCPML,new convolution perfectly matched layer,又称新卷积完全匹配层),非分裂完全匹配层吸收边界能够吸收大角度的入射波和掠射波,能够直接对二阶弹性波方程做数值模拟,减小内存使用空间,提高计算效率,更适用于实际生产需求。
附图说明
图1是本发明的基于非分裂完全匹配层吸收边界的正演模拟方法的一具体实施例流程图;
图2是本发明的具体实施例中均匀介质模型纵波源激发分裂完全匹配层(SPML)和非分裂完全匹配层(NCPML)两种不同吸收边界条件下二阶弹性波正演的水平分量波场图;(a)为SPML边界700ms时刻的波场快照;(b)为非分裂完全匹配层边界700ms时刻的波场快照;
图3是本发明的具体实施例中不同吸收边界条件下x=900m处700ms时刻水平分量波形对比图;
图4是本发明的具体实施例中层状速度模型的示意图;
图5是本发明的具体实施例中层状模型t=350ms时两种不同吸收边界条件下水平分量波场快照对比图;(a)层状模型SPML吸收边界条件下水平分量波场快照;(b)层状模型非分裂完全匹配层吸收边界条件下水平分量波场快照;
图6是本发明的具体实施例中第150道SPML吸收边界和非分裂完全匹配层吸收边界条件下的水平分量波形曲线与理论水平分量波形曲线对比图;
图7是本发明的具体实施例中图6在深度为400m至800m处局部放大水平分量波形曲线对比图;
图8是本发明的具体实施例中图6SPML吸收边界和非分裂完全匹配层吸收边界波场水平分量的数值解与理论值的误差曲线对比图;
图9是本发明的具体实施例中非分裂完全匹配层吸收边界的层状模型对应的地震记录的水平分量和垂直分量图,(a)水平分量地震记录,(b)垂直分量。
具体实施方式
为使本发明的上述和其他目的、特征和优点能更明显易懂,下文特举出较佳实施例,并配合附图所示,作详细说明如下。
图1所示为本发明的基于非分裂完全匹配层弹性波数值模拟方法的流程图。
第一步:建立二阶弹性波波动方程
通过柯西方程、几何方程以及纳维尔方程推导出二阶弹性波波动方程:
其中λ和μ为介质的拉梅常数;ρ表示介质密度;uz为波场垂直分量,ux为波场水平分量,t表示地震波传播时间,x,z分别表示水平方向和垂直方向。
CPML吸收边界x方向的复拉伸变量可以表示为:
其中,kx≥1,dx≥0,ax≥0,虚数部分dx>0表示吸收层中衰减因子,用来衰减透射波,ax,dx是新引入的衰减因子,sx表示x方向复拉伸变量;在非吸收层中,kx=1且dx=0。
采用以上复坐标变换公式(2),将二阶弹性波方程表示为:
其中和/>分别表示纵、横波速度;/>表示波场水平分量,/>为波场垂直分量;sx表示x方向复拉伸变量,sz为z方向复拉伸变量;F-1为fourier反变换。
忽略衰减因子空变影响计算出方程(3)中空间二阶偏导算子:
其中:δ(t)和H(t)分别代表狄拉克函数以及单位阶跃函数,t为时间,dx、dz为衰减透射波的衰减参数;αz、αz为滤波的频率参数;为x方向在复频域的拉伸变量,/>为x与z方向在复频域的拉伸变量,kz,kz为衰减渐消波的尺度因子。
利用迭代法代替卷积运算的思想,将时域卷积运算转换成迭代运算并计算得到时间域基于非分裂完全匹配层的二阶弹性波方程:
1)当dz/kzz≠dz/kzz时,基于非分裂完全匹配层时间域的控制方程为:
2)当dz/kzz=dz/kzz时,基于非分裂完全匹配层时间域的控制方程为:
方程(5)(6)中:是将/>卷积运算转换成迭代运算的迭代项,bz、az为转换过程中产生的两个系数,/> 类似,uz为波场垂直分量,uz为波场水平分量,vp表示纵波速度、vs表示横波速度,其余符号含义与公式(4)相同。
第三步:实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟为实现正演模拟,对连续方程进行时间二阶、空间十阶差分计算,差分方程格式如下:
时间二阶精度差分格式:
空间十阶精度差分格式:
其中a0、an均为差分系数,Δx、Δz为网格间距,u(x,z)为波场,t表示采样时间,Δt表示时间采样间隔。
离散后的弹性波方程为:
1)当dz/kzz≠dz/kzz时,基于非分裂完全匹配层时间域的控制方程的离散形式为:
2)当dz/kzz=dz/kzz时,基于非分裂完全匹配层时间域的控制方程离散形式为:
方程(10)和方程(11)中,为n+1时刻的水平分量,/>为n+1时刻垂直分量,Δt表示时间采样间隔,其它参数同公式(5)和公式(6)。
正演模拟使用的子波主频为20Hz,时间采样间隔为1ms,采样点数为1000。均匀模型震源位于网格点(100,100),层状模型震源位于(1000m,0m)处。
通过数值模拟对比非分裂完全匹配层边界和SPML边界在吸收效果、运行速度以及内存占用方面的不同,并用复杂的层状模型测试非分裂完全匹配层边界对复杂模型的适定性。首先通过层状模型对比分析非分裂完全匹配层吸收边界与SPML吸收边界的吸收效果以及计算效率;然后通过层状模型验证所推导的二阶弹性波非分裂完全匹配层吸收边界的稳定性。
图2是本发明的一具体实施例中均匀介质模型纵波源激发的水平分量波场图,模型大小为200*200,Vp=3000m/s,Vs=2000m/s,网格间距为10m。(a)为SPML边界700ms时刻的波场快照;(b)为非分裂完全匹配层边界700ms时刻的波场快照;对比图2(a)和图2(b),直观上两种吸收边界对人工边界产生的反射波都具有良好的吸收效果,在边界处没有明显的虚假反射。定量分析两种吸收边界的吸收效果,图3是图2中虚线处(在x=900m的垂线上)用两种PML吸收边界数值模拟得出的波场值与理论值的对比图,其中理论值是通过增加计算区域得到的,其计算结果等同于对边界处的透射波完全吸收。通过图2和图3可以看出与SPML吸收边界相比,非分裂完全匹配层吸收边界更接近理论值,吸收效果更好。
图4为存在多个速度界面的复杂层状速度模型。在层状模型数值模拟中,图5(a)为层状模型SPML吸收边界t=350ms时水平分量波场快照,(b)为层状模型非分裂完全匹配层吸收边界t=350ms时水平分量波场快照;图6是第150道SPML吸收边界和非分裂完全匹配层吸收边界的水平分量波形曲线与理论水平分量波形曲线对比;由图5和图6可以看出非分裂完全匹配层吸收边界与SPML吸收边界的吸收效果都较好。为进一步分析两者的区别,对图6所示曲线进行局部放大显示,并分别计算两者与理论值的误差。图7是图6在深度为400m至800m处水平分量波形曲线;图8是用图6中SPML吸收边界和非分裂完全匹配层吸收边界水平分量波场的数值解与理论值的误差曲线;通过图7和图8可以看出非分裂完全匹配层吸收边界的数值更加接近理论值,其吸收效果更好。
表1SPML边界与非分裂完全匹配层边界计算效率与占用内存数据对比
表1为均匀介质模型SPML边界与非分裂完全匹配层边界计算效率与占用内存的数据对比,通过表1的数据对比可知非分裂完全匹配层吸收边界相较于SPML吸收边界计算速度更快,所占内存更小,更适用于生产实践。图9是基于非分裂完全匹配层吸收边界的层状模型对应的地震记录,图9(a)水平分量,图9(b)垂直分量;由图9可以看出,得到的地震记录都没有明显的虚假反射波,证明了本方法对于复杂模型的数值模拟稳定性良好。

Claims (3)

1.基于非分裂完全匹配层吸收边界的正演模拟方法,其特征在于,该基于非分裂完全匹配层吸收边界的正演模拟方法包括:
步骤1,建立二阶弹性波动方程;
步骤2,利用迭代法代替卷积运算的思想构建适用于二阶弹性波方程的非分裂完全匹配层吸收边界;
步骤3,实现非分裂完全匹配层吸收边界二阶精度弹性波正演模拟;
在步骤1中,通过柯西方程、几何方程以及纳维尔方程推导出二阶弹性波波动方程:
其中λ和μ为介质的拉梅常数;ρ表示介质密度;uz为波场垂直分量,ux为波场水平分量,t表示地震波传播时间,x,z分别表示水平方向和垂直方向;
在步骤2中,采用复坐标变换,忽略衰减因子空变影响计算出二阶弹性波方程的空间二阶偏导算子:
其中:δ(t)和H(t)分别代表狄拉克函数以及单位阶跃函数,t为时间,dx、dz为衰减透射波的衰减参数;αx、αz为滤波的频率参数;为x方向在复频域的拉伸变量,/>为x与z方向在复频域的拉伸变量,kx,kz为衰减渐消波的尺度因子;
利用迭代法代替卷积运算的思想,将二阶偏导算子中的时域卷积运算转换成迭代运算计算出时间域基于非分裂完全匹配层吸收边界的二阶弹性波方程,
1)当dx/kxx≠dz/kzz时,基于非分裂完全匹配层时间域的控制方程为:
2)当dx/kxx=dz/kzz时,基于非分裂完全匹配层时间域的控制方程为:
方程(5)(6)中:是将/>卷积运算转换成迭代运算的迭代项,/>是将/>卷积运算转换成迭代运算的迭代项,/>是将/>卷积运算转换成迭代运算的迭代项,/>将/>卷积运算转换成迭代运算的迭代项,/>将/>卷积运算转换成迭代运算的迭代项,/>是将/>卷积运算转换成迭代运算的迭代项,/>将/>卷积运算转换成迭代运算的迭代项,/>将/>卷积运算转换成迭代运算的迭代项,其中上标n表示迭代次数;bx、ax为转换过程中产生的两个系数,/> uz为波场垂直分量,ux为波场水平分量,t表示地震波传播时间,vp表示纵波速度,vs表示横波速度。
2.根据权利要求1所述的基于非分裂完全匹配层吸收边界的正演模拟方法,其特征在于,在步骤3中,由连续方程的时间二阶、空间十阶有限差分计算基于非分裂完全匹配层吸收边界二阶弹性波方程的离散方程。
3.根据权利要求2所述的基于非分裂完全匹配层吸收边界的正演模拟方法,其特征在于,在步骤3中,对均匀介质模型和复杂层状模型分别开展分裂完全匹配层和非分裂完全匹配层吸收边界二阶精度弹性波方程正演模拟,通过成图验证非分裂完全匹配层吸收边界具有更好的吸收效果、更快的计算速度和对计算机内存较低的需求。
CN202010661155.5A 2020-07-10 2020-07-10 基于非分裂完全匹配层吸收边界的正演模拟方法 Active CN113917526B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010661155.5A CN113917526B (zh) 2020-07-10 2020-07-10 基于非分裂完全匹配层吸收边界的正演模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010661155.5A CN113917526B (zh) 2020-07-10 2020-07-10 基于非分裂完全匹配层吸收边界的正演模拟方法

Publications (2)

Publication Number Publication Date
CN113917526A CN113917526A (zh) 2022-01-11
CN113917526B true CN113917526B (zh) 2023-07-28

Family

ID=79231998

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010661155.5A Active CN113917526B (zh) 2020-07-10 2020-07-10 基于非分裂完全匹配层吸收边界的正演模拟方法

Country Status (1)

Country Link
CN (1) CN113917526B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107479092A (zh) * 2017-08-17 2017-12-15 电子科技大学 一种基于方向导数的频率域高阶声波方程正演模拟方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7856329B2 (en) * 2006-04-10 2010-12-21 Tim Niebauer Method and apparatus for processing an under-sampled chirped sinusoidal waveform using a complex-heterodyne
CN102253415B (zh) * 2011-04-19 2013-03-20 中国石油大学(华东) 基于裂缝等效介质模型的地震响应模式建立方法
CN104459791A (zh) * 2014-12-15 2015-03-25 中国石油大学(华东) 一种基于波动方程的小尺度大模型正演模拟方法
US10520618B2 (en) * 2015-02-04 2019-12-31 ExxohnMobil Upstream Research Company Poynting vector minimal reflection boundary conditions
US10317546B2 (en) * 2015-02-13 2019-06-11 Exxonmobil Upstream Research Company Efficient and stable absorbing boundary condition in finite-difference calculations
CN106777472B (zh) * 2016-11-16 2020-05-22 西安理工大学 基于拉盖尔多项式的减少分裂误差的完全匹配层实现方法
CN107944214B (zh) * 2017-11-27 2020-11-10 河北工业大学 笛卡尔坐标系下各向异性完全匹配层截断边界的实现方法
CN111208563B (zh) * 2020-02-18 2021-08-06 吉林大学 一种非分裂完全匹配层吸收边界方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107479092A (zh) * 2017-08-17 2017-12-15 电子科技大学 一种基于方向导数的频率域高阶声波方程正演模拟方法

Also Published As

Publication number Publication date
CN113917526A (zh) 2022-01-11

Similar Documents

Publication Publication Date Title
CN103293552B (zh) 一种叠前地震资料的反演方法及系统
Carcione Theory and modeling of constant-Q P-and S-waves using fractional time derivatives
Zhu et al. Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians
KR101948509B1 (ko) 지구 물리학적 데이터의 반복 반전의 아티팩트 감소
US20140372043A1 (en) Full Waveform Inversion Using Perfectly Reflectionless Subgridding
CN110456417B (zh) 一种地震数据多次波压制方法
CN112327358B (zh) 一种粘滞性介质中声波地震数据正演模拟方法
Chen et al. A matrix-transform numerical solver for fractional Laplacian viscoacoustic wave equation
CN109946742B (zh) 一种TTI介质中纯qP波地震数据模拟方法
CN110018517A (zh) 一种多尺度地面微地震逆时干涉定位方法
Ozgun Konca et al. Kinematic inversion of physically plausible earthquake source models obtained from dynamic rupture simulations
CN110967749A (zh) 一种vsp地震资料频变q值估计与反q滤波方法
Qin et al. The implementation of an improved NPML absorbing boundary condition in elastic wave modeling
AU2017366757A1 (en) Seismic acquisition geometry evaluation using full-waveform inversion
CN109738950A (zh) 基于三维稀疏聚焦域反演的噪声型数据一次波反演方法
CN111665556B (zh) 地层声波传播速度模型构建方法
CN110658558A (zh) 吸收衰减介质叠前深度逆时偏移成像方法及系统
CN109164488B (zh) 一种梯形网格有限差分地震波场模拟方法
AU2015290248A1 (en) System and method for rock property estimation of subsurface geologic volumes
CN113655522A (zh) 频率域地震弱信号的增强方法
CN113917526B (zh) 基于非分裂完全匹配层吸收边界的正演模拟方法
Wang et al. Rock Fracture Monitoring Based on High‐Precision Microseismic Event Location Using 3D Multiscale Waveform Inversion
CN113866823A (zh) 一种粘声各向异性介质中的正演成像方法
CN116882214B (zh) 基于dfl粘弹性方程的瑞雷波数值模拟方法及系统
CN114114406B (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