CN106680871B - 一种基于p波到时与波形振幅的震源定位方法 - Google Patents

一种基于p波到时与波形振幅的震源定位方法 Download PDF

Info

Publication number
CN106680871B
CN106680871B CN201611099983.4A CN201611099983A CN106680871B CN 106680871 B CN106680871 B CN 106680871B CN 201611099983 A CN201611099983 A CN 201611099983A CN 106680871 B CN106680871 B CN 106680871B
Authority
CN
China
Prior art keywords
mrow
msub
msup
wave
amplitude
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
CN201611099983.4A
Other languages
English (en)
Other versions
CN106680871A (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
Original Assignee
Central South University
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 filed Critical Central South University
Priority to CN201611099983.4A priority Critical patent/CN106680871B/zh
Publication of CN106680871A publication Critical patent/CN106680871A/zh
Application granted granted Critical
Publication of CN106680871B publication Critical patent/CN106680871B/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/30Analysis
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/65Source localisation, e.g. faults, hypocenters or reservoirs

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)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种基于P波到时与波形振幅的震源定位方法,包括以下步骤:步骤1:通过微震监测系统获取震源参数:步骤2:在未知P波波速v的情况下,求解一次震源坐标(x0,y0,z0)及相应速度v;步骤3:计算波形振幅趋势线的斜率步骤4:以logistic函数拟合波形振幅趋势线的斜率K与每个震源的P波波速v之间的函数关系式;步骤5:计算修正后的各个传感器对应的P波波速步骤6:利用修正后的各个传感器对应的P波波速求解二次震源坐标(x,y,z)。本发明方法具有适用性强、准确性高等特点。

Description

一种基于P波到时与波形振幅的震源定位方法
技术领域
本发明涉及一种基于P波到时与波形振幅的震源定位方法,结果明确,适用性强,准确率高。
背景技术
微震监测作为一种地下工程安全监测的重要技术手段之一,已经广泛被应用于矿山安全、水电工程等领域。而微震震源空间位置的确定以及微震震源定位精度的提高,一直是微震监测技术的重要内容。但是,微震与爆破事件波形较为复杂,S波的首波通常夹杂在P波的尾波中,使得S波到时不仅难以辨识,且易受操作人员专业知识和经验的影响,极大降低了S波拾取准确度及数量,严重制约微震监测实时分析的效果。
震源定位方法到目前为止有十余种之多,主要包括几何方法、物理方法和数学方法等。Geiger定位方法是通过迭代计算使时间残值最小化。在此基础上,Buland对其进行修改并发展出基于经典的最小二乘法定位计算,使其能更好适应地震定位。Tatantola提出了Bayesian定位方法的严格公式和解。Spence针对时差定位技术模型的不确定导致走时异常进行了有意义的研究。Crosson应用联合测定方法(SSH)对局部事件进行了定位研究。以上定位方法,都是以预先测定平均速度或给出平均速度模型为前提,故在现场进行定位时,平均速度测量的准确性直接影响着定位精度。而且,对S波到时的自动准确拾取仍存在难题,需要依靠具有经验的专业工作人员人工拾取,定位准确性与精度存在较大的随机性。
可见现有的微震与爆破事件震源定位方法存在较大的局限,需要研究一种适用性强、准确率高的定位方法。
发明内容
本发明所要解决的技术问题是提供一种基于P波到时与波形振幅的震源定位方法,该定位方法适用性强、准确率高。
发明的技术解决方案如下:
一种基于P波到时与波形振幅的震源定位方法,包括以下步骤:
步骤1:通过微震监测系统获取震源参数;
震源参数包括:各传感器坐标(xi,yi,zi),其中i为传感器编号,i=1,2,…,n;震源起振时刻t0,各传感器的P波到时tPi,各传感器的波形振幅趋势线的最大振幅MaxAmpi及其对应时刻tMi
步骤2:根据步骤1获取的初步求解震源坐标(x,y,z)及P波波速v,得到一次定位的震源坐标(x0,y0,z0)和初步求解的P波波速v0
步骤3:计算各传感器i(i=1,2,…,n)的波形振幅趋势线的斜率Ki
步骤4:基于步骤2得到的v0和步骤3得到所有传感器的波形振幅趋势线的斜率Ki(i=1,2,…,n),通过logistic函数拟合得到传感器的波形振幅趋势线的斜率K与震源的P波波速vc之间的函数关系式K(vc):
其中,r为传感器的波形振幅趋势线的斜率随速度的增长速率(r是logistic函数的一个参数,经数值分析软件拟合数据之后可直接得出),Kmax为所有传感器的波形振幅趋势线的斜率中的最大值,即Kmax=max{Ki(i=1,2,…,n)};K0为所有传感器的波形振幅趋势线的斜率中的最小值,即Kmax=min{Ki(i=1,2,…,n)};
步骤5:计算修正后的传感器i对应的P波波速
保持传感器i的波形振幅趋势线的斜率Ki不变,将Ki作为K(vc)代入公式二,计算相应的vc值,即为修正后的传感器i对应的P波波速
步骤6:二次定位震源坐标(x,y,z);
将修正后的传感器i对应的P波波速代入以下目标函数:
其中,min表示求最小值;Δtij=tPi-tPj
求解目标函数,其最优解即为二次定位的震源坐标(x,y,z)。
所述步骤2具体为:
建立以下目标函数:
其中,min表示求最小值;Δtij=tPi-tPj
由此求得使得全部理论值与实测值Δtij的偏离平方和最小时的震源坐标(x,y,z)及P波波速v,即一次定位的震源坐标(x0,y0,z0)和初步求解的P波波速v0
有益效果:
本发明选取震源振动波形图中的最易准确拾取的参数:震源起振时刻t0,各传感器P波到时tPi,最大振幅MaxAmpi及其对应时刻tMi,优化定位所用的P波波速,推导求差式非线性拟合函数用于震源定位,事实证明其结果明确,适用性强,准确率高。
本发明基于P波到时与波形振幅对震源坐标进行确定,减轻了人工拾取S波的数据处理压力,避免了因S波拾取不准确而带来的误差,提高了定位精度与准确性。此方法具有适用性强、准确性高等特点。
附图说明
图1为本发明流程图。
具体实施方式
下面将结合实际算例,对本发明提出的一种基于P波到时与波形振幅的震源定位方法作进一步说明。本发明的思想描述如下:本发明通过微震监测系统获取主要震源参数;利用速度—距离公式和两点间距离公式得到震源与任意传感器i、j之间的等量关系,求得一次震源坐标;再通过优化P波波速提高准确性,推导求差式非线性拟合函数用于震源定位,计算出震源的准确坐标。此方法具有适用性强、准确率高的优点。
步骤1:通过微震监测系统获取震源参数;
震源参数包括:各传感器坐标(xi,yi,zi),其中i为传感器编号,i=1,2,…,n;震源起振时刻t0,各传感器的P波到时tPi,各传感器的波形振幅趋势线的最大振幅MaxAmpi及其对应时刻tMi
步骤2:根据步骤1获取一次定位的震源坐标(x,y,z)及初步求解的P波波速v,分别记为(x0,y0,z0)和v0;;
设li为传感器i至震源的距离,由两点间距离公式可得:
由速度—距离公式可知,传感器i的P波到时的理论值ti为:
即:
由式(4)可确定任意2个不同的传感器i和j的P波到时之差的理论值为:
计算传感器i和j的P波到时之差的实测值Δtij
Δtij=tPi-tPj (5)
全部理论值与实测Δtij的偏离平方和可描述全部理论值与实测值的偏离程度,其偏离越小,则相应的一次定位的震源坐标(x0,y0,z0)及初步求解的P波波速v0与其真实值的拟合度越好,准确度越高,因此,建立以下目标函数:
其中,min表示求最小值,
由此求得使得全部理论值与实测值Δtij的偏离平方和最小时的震源坐标(x,y,z)及P波波速v,即一次定位的震源坐标(x0,y0,z0)和初步求解的P波波速v0;;
步骤3:计算传感器i(i=1,2,…,n)的波形振幅趋势线的斜率Ki
步骤4:基于步骤2初步求解的P波波速v0和步骤3得到所有传感器的波形振幅趋势线的斜率Ki(i=1,2,…,n),通过logistic函数拟合得到传感器的波形振幅趋势线的斜率K与震源的P波波速vc之间的函数关系式K(vc):
其中,r为传感器的波形振幅趋势线的斜率随速度的增长速率(r是logistic函数的一个参数,经数值分析软件拟合数据之后可直接得出),Kmax为所有传感器的波形振幅趋势线的斜率中的最大值,即Kmax=max{Ki(i=1,2,…,n)};K0为所有传感器的波形振幅趋势线的斜率中的最小值,即Kmax=min{Ki(i=1,2,…,n)};
步骤5:计算修正后的传感器i对应的P波波速
保持传感器i的波形振幅趋势线的斜率Ki不变,将Ki作为K(vc)代入函数关系式(8),计算相应的vc值,即为修正后的传感器i对应的P波波速
步骤6:再次求解震源坐标(x,y,z);
将修正后的传感器i对应的P波波速作为v,代入目标函数(6)有:
上式必定取得非负值,因此其最小值一定存在。求解(10),其最优解即为二次定位的震源坐标(x,y,z)。
本发明选取震源振动波形图中的最易准确拾取的参数:震源起振时刻t0,各传感器P波到时tPi,最大振幅MaxAmpi及其对应时刻tMi,优化定位所用的P波波速,推导求差式非线性拟合函数用于震源定位,事实证明其结果明确,适用性强,准确率高。实施实例:
由步骤2可得一次定位震源坐标如下所示:
表1 5组震源的实际坐标和一次定位坐标
由步骤3、步骤4、步骤5可得表2如下所示:
表2 5组震源的K-v函数关系式
由步骤6可求得二次震源坐标如下表所示
表3 5组震源的实际坐标和二次定位坐标
由表3可知,经速度修正后,二次定位坐标精度明显高于一次定位坐标精度。
以上所述仅为本发明的实施例而已,并不用以限制本发明,凡在本发明精神和原则之内,所作任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于P波到时与波形振幅的震源定位方法,其特征在于,包括以下步骤:
步骤1:通过微震监测系统获取震源参数;
震源参数包括:各传感器坐标(xi,yi,zi),其中i为传感器编号,i=1,2,…,n;震源起振时刻t0,各传感器的P波到时tPi,各传感器的波形振幅趋势线的最大振幅MaxAmpi及其对应时刻tMi
步骤2:根据步骤1获取的震源参数初步求解震源坐标(x,y,z)及P波波速v,得到一次定位的震源坐标(x0,y0,z0)和初步求解的P波波速v0
步骤3:计算各传感器i的波形振幅趋势线的斜率Ki
步骤4:基于步骤2初步求解的P波波速v0和步骤3得到所有传感器的波形振幅趋势线的斜率Ki,通过logistic函数拟合得到传感器的波形振幅趋势线的斜率K与震源的P波波速vc之间的函数关系式K(vc):
其中,r为传感器的波形振幅趋势线的斜率随速度的增长速率,Kmax为所有传感器的波形振幅趋势线的斜率中的最大值,K0为所有传感器的波形振幅趋势线的斜率中的最小值;
步骤5:计算修正后的传感器i对应的P波波速
将Ki作为K(vc)代入公式二,计算相应的vc值,即为修正后的传感器i对应的P波波速
步骤6:二次定位震源坐标(x,y,z);
将修正后的传感器i对应的P波波速代入以下目标函数:
<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>min</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;Delta;t</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mfrac> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>x</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>y</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>z</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <msubsup> <mi>v</mi> <mi>i</mi> <mi>C</mi> </msubsup> </mfrac> <mo>+</mo> <mfrac> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>x</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>z</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <msubsup> <mi>v</mi> <mi>i</mi> <mi>C</mi> </msubsup> </mfrac> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow>
其中,min表示求最小值;Δtij=tPi-tPj
求解目标函数,其最优解即为二次定位的震源坐标(x,y,z)。
2.根据权利要求1所述的基于P波到时与波形振幅的震源定位方法,其特征在于,所述步骤2具体为:
建立以下目标函数:
<mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>min</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>&amp;Delta;t</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mi>&amp;Delta;</mi> <msub> <mover> <mi>t</mi> <mo>^</mo> </mover> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow>
其中,min表示求最小值;Δtij=tPi-tPj
<mrow> <mi>&amp;Delta;</mi> <msub> <mover> <mi>t</mi> <mo>^</mo> </mover> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mfrac> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <mi>x</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>y</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>z</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mi>v</mi> </mfrac> <mo>-</mo> <mfrac> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>x</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>y</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>j</mi> </msub> <mo>-</mo> <mi>z</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mi>v</mi> </mfrac> <mo>;</mo> </mrow>
由此求得使得全部理论值与实测值Δtij的偏离平方和最小时的震源坐标(x,y,z)及P波波速v,即一次定位的震源坐标(x0,y0,z0)和初步求解的P波波速v0
3.根据权利要求1所述的基于P波到时与波形振幅的震源定位方法,其特征在于,所述公式二中的r经数值分析软件拟合数据之后可直接得出。
CN201611099983.4A 2016-12-05 2016-12-05 一种基于p波到时与波形振幅的震源定位方法 Active CN106680871B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201611099983.4A CN106680871B (zh) 2016-12-05 2016-12-05 一种基于p波到时与波形振幅的震源定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201611099983.4A CN106680871B (zh) 2016-12-05 2016-12-05 一种基于p波到时与波形振幅的震源定位方法

Publications (2)

Publication Number Publication Date
CN106680871A CN106680871A (zh) 2017-05-17
CN106680871B true CN106680871B (zh) 2017-09-26

Family

ID=58867457

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201611099983.4A Active CN106680871B (zh) 2016-12-05 2016-12-05 一种基于p波到时与波形振幅的震源定位方法

Country Status (1)

Country Link
CN (1) CN106680871B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107144876A (zh) * 2017-05-24 2017-09-08 辽宁大学 一种基于时窗的煤矿冲击地压震级及震源定位方法
CN109597125B (zh) * 2018-11-27 2020-10-16 湖北海震科创技术有限公司 一种基于p波到时与最大振幅波形的微震源定位方法
CN109521221B (zh) * 2018-12-10 2020-09-29 东北大学 一种钻爆法施工硬岩隧道微震波波速实时获取方法
CN109828236A (zh) * 2019-02-14 2019-05-31 中南大学 一种含空区复杂结构中的微震/声发射源定位方法
CN111457252B (zh) * 2020-06-01 2022-07-19 安徽理工大学 一种基于振动波的燃气管道泄漏定位方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102778668A (zh) * 2012-07-23 2012-11-14 中煤科工集团西安研究院 矿山被动震源快速精确定位方法
CN103389489B (zh) * 2013-07-31 2015-03-04 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于大斜度井的微地震监测定位方法
CN105759311A (zh) * 2016-01-25 2016-07-13 西南交通大学 一种近实时地震震源位置定位方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3463677B2 (ja) * 2001-10-04 2003-11-05 独立行政法人防災科学技術研究所 震源位置の決定法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102778668A (zh) * 2012-07-23 2012-11-14 中煤科工集团西安研究院 矿山被动震源快速精确定位方法
CN103389489B (zh) * 2013-07-31 2015-03-04 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于大斜度井的微地震监测定位方法
CN105759311A (zh) * 2016-01-25 2016-07-13 西南交通大学 一种近实时地震震源位置定位方法

Also Published As

Publication number Publication date
CN106680871A (zh) 2017-05-17

Similar Documents

Publication Publication Date Title
CN106680871B (zh) 一种基于p波到时与波形振幅的震源定位方法
Madhukar et al. State estimation using extended kalman filter and unscented kalman filter
CN103900574B (zh) 一种基于迭代容积卡尔曼滤波姿态估计方法
CN105159084B (zh) 一种带干扰观测器的机械手神经网络控制系统及方法
CN104697436B (zh) 一种基于傅里叶级数的圆感应同步器误差模型分析方法
CN104964685B (zh) 一种手机运动姿态的判定方法
CN104197927A (zh) 水下结构检测机器人实时导航系统及方法
WO2022134144A1 (zh) 一种机器人质心规划方法、装置、可读存储介质及机器人
CN103604430A (zh) 一种基于边缘化ckf重力辅助导航的方法
CN106323334A (zh) 一种基于粒子群优化的磁力计校准方法
CN104865597A (zh) 一种深度域层速度初始模型的建模方法
CN102540252A (zh) 基于互相关的高精度中值叠加方法
CN109341725A (zh) 行星接近段导航性能快速评估方法
US20130079949A1 (en) Inclination detection systems and methods
CN108303095A (zh) 适用于非高斯系统的鲁棒容积目标协同定位方法
CN106091958A (zh) 基于弓高弦长法测量圆弧工件半径的方法
CN103438872A (zh) 一种基于大坝三维前方交会测量的内外业一体化系统
CN105093280A (zh) 表层模型对地震数据影响的低频与高频成分的分解方法
CN107644119B (zh) 一种基于三维扫描点云的半孔率自动计算方法
CN104765476A (zh) 手写轨迹生成方法和装置
WO2023202157A1 (zh) 铲斗坐标标定方法和装置、更新方法和设备、挖掘机
CN103745073A (zh) 一种边坡三维变形预测方法
CN109866217A (zh) 机器人里程定位方法、装置、终端设备及计算机存储介质
CN106092141B (zh) 一种改善相对位置传感器性能的方法及装置
CN108919313A (zh) 利用最优数值导数的gnss多普勒观测值生成方法

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