CN105426339B - 一种基于无网格法的线源时域电磁响应数值计算方法 - Google Patents
一种基于无网格法的线源时域电磁响应数值计算方法 Download PDFInfo
- Publication number
- CN105426339B CN105426339B CN201510751170.8A CN201510751170A CN105426339B CN 105426339 B CN105426339 B CN 105426339B CN 201510751170 A CN201510751170 A CN 201510751170A CN 105426339 B CN105426339 B CN 105426339B
- Authority
- CN
- China
- Prior art keywords
- line source
- node
- domain
- equation
- electromagnetic
- 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
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
- G06F17/147—Discrete orthonormal transforms, e.g. discrete cosine transform, discrete sine transform, and variations therefrom, e.g. modified discrete cosine transform, integer transforms approximating the discrete cosine transform
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Discrete Mathematics (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种基于无网格法的线源时域电磁响应数值计算方法,尤其是可以克服传统数值计算方法中对于网格的依赖,适用于复杂地形下时域电磁探测的数值模拟。本发明基于瞬变电磁法满足的控制方程和定解条件,建立了二维线源边值问题的泛函,利用罚因子法加载本质边界条件,提出旁轴近似方程消除截断边界处的反射波,采用Crack‑Nicolson格式进行时间离散,得到递推方程。利用等参单元思想将局部坐标中形状规则的单元离散为节点任意分布的不规则求解对象。采用LU分解方法求解递推方程,最终得到求解区域内各个节点的场值。计算结果表明,该方法形函数光滑性好,模拟精度高,最大误差不超过1×10‑3,实现了电磁法高精度的数值计算。
Description
技术领域
本发明涉及一种地球物理勘探领域的电磁场数值计算方法,尤其是基于无网格法的线源时域电磁响应数值计算方法。
背景技术
时间域电磁法(Time domain electromagnetic methods)或称瞬变电磁法(Transient electromagnetic methods),是一种建立在电磁感应原理基础上的人工源电磁探测方法。目前该方法被广泛应用于金属矿勘探、煤矿水文地质调查和工程勘查等领域。
电磁场数值模拟是数据处理和反演解释的基础,在地球物理数值计算中具有基础性、全局性作用。现有瞬变电磁数值计算方法主要有积分方程法(IEM),有限差分法(FDM),有限单元法(FEM)等。
CN201410125866.5公开了一种频率域正演方法及装置,该方法通过建立17点格式的差分公式,构造稀疏矩阵,读入子波参数和速度模型,频率循环得到单频波场,对所有频率波场求反傅里叶变换,得到正演结果。该方法增加了差分公式的适用性,通过加载PML(最佳匹配层)边界条件,防止边界反射的波场回折到地表。
CN201110459223.0公开了一种电磁场仿真分析方法,该方法采用棱边单元对三维模型进行网格剖分,引入满足inf-sub条件的标量乘子空间,加入有约束条件的集成后总体矩阵方程求解,得到电磁场分析结果。
美国专利US11756384公开了一种基于三维时域正演建模基础上的反演方法。该方法基于梯度的波形反演理论,通过匹配建模的数据来估计模型参数,在时域建模之后通过离散傅里叶变换(FFT)转换到频率域,完成了三维频域波形反演。
以上所述方法均是基于网格实现的,在计算的过程中都不可避免地涉及到网格的剖分问题,对于单元的形状也有一定的要求,特别是对于物性参数分布复杂和几何特征分布不规则的地电模型适应性差。
发明内容
本发明所要解决的技术问题在于提供一种基于无网格法的线源时域电磁响应数值计算方法,可以克服传统正演方法中对于网格的依赖,适用于复杂地形下时域电磁探测的数值模拟。
本发明是这样实现的,一种基于无网格法的线源时域电磁响应数值计算方法包括:
1)在计算区域内进行节点、电性参数、支持域、背景网格、初始场值设置;
2)在计算区域内对背景网格进行等参数变换,在每个背景单元中分别计算高斯积分点及其权系数;
3)对所有背景网格进行外循环,对所有高斯积分点内循环,搜索背景网格中局部定义域内的有效节点并计算局部定义域内节点处形函数;
4)加载强加边界条件得到线性方程组,加载吸收边界条件得到边界阻尼矩阵,将边界阻尼矩阵加入到线性方程组中;
5)采用LU分解法求解线性方程组,得到各个节点各个时刻的场值。
进一步地,步骤3中局部定义域内利用滑动最小二乘法建立形函数。
进一步地,步骤4中利用罚因子法处理计算模型的强加边界条件,建立和瞬变电磁边值问题的等价泛函数,求取泛函数的变分,得到刚度矩阵K,阻尼矩阵K'和右端项S,形成线性方程组:
进一步地,步骤4中根据旁轴近似方程加载吸收边界条件,将左行波作用于右行入射的波,右行波作用于左行入射的波,下行波作用于上行入射的波,相互抵消后消除截断边界处的反射波,得到边界阻尼矩阵:
N1和N2均为无网格法边界处两个节点的形函数,nx和ny为边界外法线的方向余弦,v为地下介质中的电磁波传播速度。
进一步地,步骤5中线性方程组的时间变量的离散采用Crack-Nicolson格式,令K′+F=K”,最终形成的递推方程为:
(2K”+ΔtKt+Δt)Et+Δt=ΔtSt+Δt+StΔt+2K”Et-ΔtKtEt
其中Et为当前时刻的场值,Et+Δt为待求的下一时刻的场值,S为线性方程组的右端项,Δt为离散的时间步长,K为刚度矩阵。
本发明与现有技术相比,有益效果在于:本发明针对传统的网格化数值计算方式的不足,提出了无网格法模拟的方案。以二维线源瞬变电磁法为例,利用变分原理详细推导了等价线性方程组的形成过程,采用Crack-Nicolson方案离散时间变量,形成最终的递推方程,求取了线源各个节点各个时刻的电磁响应,完整地展现了无网格法电磁数值模拟的全部过程。由于无网格法摆脱了传统方法中对于网格的依赖,理论上节点可以任意布置,为含有起伏地形等复杂情况下电磁法的数值模拟提供新的参考。在实践中同时发现无网格法形函数光滑性好,模拟精度高,为实现电磁法高精度反演打下了坚实的基础。
附图说明
图1是无网格法高斯点、场节点、影响域,局部定义域与背景单元示意图;
图2是二维线源时域电磁探测响应计算模拟示意图;
图3是截断边界处的旁轴近似方程示意图;
图4是任意四边形下的等参数变换示意图;
图5是二维线源时域电磁响应计算方法流程图;
图6是本发明一个实施例节点处的无网格法数值解与理论解对比图;
图7是本发明一个实施例一条剖面中无网格法数值解与理论解对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
实施例
参见图5结合图1所示,一种基于无网格法的线源时域电磁响应数值计算方法,包括:
1)设置模型参数,包括节点设置、电性参数设置、背景网格设置、初始场值设置;
如图1所示的无网格法高斯点、场节点、影响域ΩI,局部定义域Ωx与背景单元示意图,包括背景积分网络1中的求解域Ω,在一个背景积分网络1的求解域Ω内包括用黑点表示的节点3,以及用空心圆表示的高斯积分点或计算点2,影响域ΩI是指一个节点的影响区域,影响域ΩI的中心为场节点,而局部定义域Ωx则为计算点进行最小二乘拟合所选择的区域,中心常常为一个计算点。在实际计算中常用节点的影响域取代计算点的局部定义域,即在构造一个计算点的无网格形函数时,若计算点位于场节点的影响域中,则该场节点属于构造该形函数的场节点,也即如果场节点的影响域包含该计算点,则该场节点参与构造该计算点的形函数。
2)在计算区域内对背景网格进行等参数变换,在每个背景单元中分别计算高斯积分点及其权系数;
3)对所有背景网格进行循环;
4)对该背景网格的所有高斯积分点循环;
判断高斯积分点是否位于求解域Ω外,如是,则忽略该高斯积分点,处理下个高斯点;
如高斯积分点位于求解域Ω内,计算此计算点的形函数,带入求积公式中计算积分,并将其组装到线性方程组的刚度矩阵和阻尼矩阵中;
5)搜索局部定义域内的有效节点;
6)计算局部定义域内节点处形函数,加载强加边界条件得到线性方程组,加载吸收边界条件得到边界阻尼矩阵,将边界阻尼矩阵加入到线性方程组中,利用Crack-Nicolson方法离散时间变量,得到最终的递推方程;
7)采用LU分解法求解递推方程,得到各个节点各个时刻的场值。
步骤6中,局部支持域Ωx内滑动最小二乘法建立形函数:
如图1所示,在计算点的局部定义域Ωx内,构造场函数的近似形式
pT(x)代表基函数,m为基函数的维数,a(x)为m维系数矢量,基函数pT(x)通常为Pascal三角形所决定的单项式以保证其最小完备性。
在局部范围内构造加权离散L2范数J得:
ωi(x)=w(x-xi)称为权函数,n为权函数非零域内的节点个数。以J取最小值得关于a(x)的支配方程:
展开可得如下线性方程:
A(x)a(x)=B(x)U
其中:
将a(x)代入场函数的近似形式(1),可得:
Ni(x)即为节点i处的形函数,它是坐标(x,y)的函数。写成矩阵形式即为:
N(x)=PT(x)A-1(x)B(x) (3)
步骤6中,二维线源瞬变电磁无网格法离散:
如图2所示,在二维空间中建立x-z坐标系,z轴以向下为正,二维线源瞬变电磁法的边值问题可以如下描述:
E(x,z,t)为待求的场量,Γ1为上边界,t0为初始时刻。μ是地下媒质的磁导率,σ是媒质的电导率,s(t)为上边界的场值。利用罚因子法处理上边界(强加边界条件),建立和上述边值问题的等价泛函:
a为罚因子,为一固定的常数。将E(x,z,t)=N(x,z)E(t)代入(4)中得:
其中,N(x,z)指的是形函数。
将最后一项展开,可得:
求取泛函I对ET(t)的偏导数,可得:
令可得:
可以简化为如下线性方程组:
其中:
其中,K称为刚度矩阵,K'为阻尼矩阵,S为右端项。
步骤6中,吸收边界条件的加载:
如图3所示,在均匀各向同性介质中,根据StanfordUniversity Claerbout(1985)推导出来的旁轴近似方程,下行波可以表示为:
左行波方程为:
右行波方程为:
由法向导数的定义知:
其中,nx和ny为边界外法线的方向余弦,v为地下介质中的电磁波传播速度。
将左行波作用于右行入射的波,右行波作用于左行入射的波,下行波作用于上行入射的波,相互抵消后实现了消除截断边界处的反射波。利用高斯定理可得:
其中,Γl、Γd、Γr分别表示左边界、底边界和右边界。这里没有上边界是因为上边界为自由边界条件。N是指形函数。
将左行波,右行波和下行波方程分别代入到(6)式,可以得到:
则边界阻尼矩阵为:
将边界阻尼矩阵F加入到线性方程组(5)中,便得到了完整的二维线源瞬变电磁无网格法数值模拟的等价方程:
在时域内进行时间离散:
令K′+F=K”,当时,称为C-N(Crank-Nicolson)格式,最终的递推方程为:
(2K”+ΔtKt+Δt)Et+Δt=ΔtSt+Δt+StΔt+2K”Et-ΔtKtEt (8)
其中Et为当前时刻的场值,Et+Δt为待求的下一时刻的场值,S为线性方程组的右端项,Δt为离散的时间步长,K为刚度矩阵。
初始场值由边值条件给定,由方程(8)实现递推从而得到各个时刻各个节点处的场值。
步骤2中,任意背景网格积分的等参数变换:
对于任意设置的背景网格积分区域如任意四边形,可以利用等参单元方法确定高斯积分点和权重。
如图4所示,对于四边形的积分区域,建立(ξ,η)→(x,y)的坐标映射,得到一阶偏导数的雅克比矩阵为:
其中,
对于二维单元积分的高斯积分有:
式中,Hi,Hj,(ξj,ηi)分别为一维情况下的高斯积分权系数和积分点,|J(ξj,ηi)|为对应的雅克比行列式值。
如图2所示的实施例,无限长导线源置于点(0,0)处,线源电流方向平行于地质体走向,影响域无量纲尺寸α=1,权函数采用三次样条函数,节点横向不均匀分布,源附近较密,最小间距1m,纵向节点均匀分布,间距为10m,均匀半空间电阻率为100Ω·m。
图6为在计算区域内抽取的任意一点的数值解与理论解的对比图,图7为选取的地下任一深度内的剖面曲线,两者高度吻合,最大误差不超过1×10-3,充分验证了无网格法高精度模拟的特点,为电磁法的数值模拟提供了新的思路和方法。
以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。
Claims (1)
1.一种基于无网格法的线源时域电磁响应数值计算方法,其特征在于:
包括如下的步骤:
1)在计算区域内进行节点、电性参数、支持域、背景网格、初始场值设置;
2)在计算区域内对背景网格进行等参数变换,在每个背景单元中分别计算高斯积分点及其权系数;
3)对所有背景网格进行外循环,对所有高斯积分点内循环,搜索背景网格中局部定义域内的有效节点并计算局部定义域内节点处形函数;
4)从线源瞬变电磁满足的控制方程和边界条件出发,
E(x,z,t)为待求的场量,Γ1为上边界,t0为初始时刻,μ是地下媒质的磁导率,σ是媒质的电导率,s(t)为上边界的场值;
建立和时域电磁边值问题的等价泛函数,利用罚因子法处理计算模型的强加边界条件,求取泛函数的变分,得到刚度矩阵K,阻尼矩阵K'和右端项S,E为待求的场值,形成线性方程组:
其中:
其中,N指的是形函数矩阵,a为罚因子;
吸收边界条件采用旁轴近似方程加载,得到边界阻尼矩阵,
N1和N2均为无网格法边界处两个节点的形函数,nx和ny为边界外法线的方向余弦,v为地下介质中的电磁波传播速度;
将边界阻尼矩阵加入到线性方程组中;得到了无网格法线源时域电磁响应数值计算的等价方程组:
5)线性方程组的时间变量的离散采用Crack-Nicolson格式,令K′+F=K”,最终形成的递推方程为:
(2K”+ΔtKt+Δt)Et+Δt=ΔtSt+Δt+StΔt+2K”Et-ΔtKtEt
其中Et为当前时刻的场值,Et+Δt为待求的下一时刻的场值,S为线性方程组的右端项,Δt为离散的时间步长,K为刚度矩阵;
采用LU分解法求解线性方程组,得到各个节点各个时刻的场值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510751170.8A CN105426339B (zh) | 2015-11-06 | 2015-11-06 | 一种基于无网格法的线源时域电磁响应数值计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510751170.8A CN105426339B (zh) | 2015-11-06 | 2015-11-06 | 一种基于无网格法的线源时域电磁响应数值计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105426339A CN105426339A (zh) | 2016-03-23 |
CN105426339B true CN105426339B (zh) | 2018-05-29 |
Family
ID=55504554
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510751170.8A Active CN105426339B (zh) | 2015-11-06 | 2015-11-06 | 一种基于无网格法的线源时域电磁响应数值计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105426339B (zh) |
Families Citing this family (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105808968B (zh) * | 2016-04-13 | 2018-07-06 | 吉林大学 | 一种时域航空电磁数值模拟中c-pml边界条件加载方法 |
CN108073732A (zh) * | 2016-11-10 | 2018-05-25 | 中国石油化工股份有限公司 | 获得稳定的近最佳匹配层吸收边界条件的方法 |
CN107798190B (zh) * | 2017-10-26 | 2021-04-09 | 吉林大学 | 复杂地形下的时域地空瞬变电磁三维数值模拟方法 |
CN107991711B (zh) * | 2017-11-27 | 2019-04-19 | 吉林大学 | 航空时域电磁三维条状随机断裂带模型建立及判别方法 |
CN108108578A (zh) * | 2018-01-30 | 2018-06-01 | 南京理工大学 | 基于无网格法的fg-grc板屈曲载荷因子的数值算法 |
CN108873084B (zh) * | 2018-05-10 | 2019-10-08 | 中南大学 | 一种基于单位分解积分的直流电阻率无单元正演方法 |
CN108875173B (zh) * | 2018-06-05 | 2022-04-26 | 哈尔滨工业大学深圳研究生院 | 一种无网格伽辽金支持域节点选取方法 |
CN109117458A (zh) * | 2018-06-28 | 2019-01-01 | 浙江省电力有限公司电力科学研究院 | 一种基于改进无网格法的直流设备积污特性计算方法 |
CN112649859B (zh) * | 2019-10-12 | 2024-03-22 | 中国石油化工股份有限公司 | 一种地震波速度自适应无网格场节点建立方法及系统 |
CN112287544B (zh) * | 2020-10-29 | 2022-02-25 | 吉林大学 | 一种无网格法的二维分形目标体频域电磁数值模拟方法 |
CN112597649B (zh) * | 2020-12-22 | 2022-04-29 | 国网湖北省电力有限公司电力科学研究院 | 一种强弱耦合的无网格静电场数值计算方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101944144A (zh) * | 2010-08-30 | 2011-01-12 | 陈玉君 | 一种基于无网格的布类仿真方法 |
CN102262699A (zh) * | 2011-07-27 | 2011-11-30 | 华北水利水电学院 | 基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法 |
CN103970717A (zh) * | 2014-05-08 | 2014-08-06 | 中国人民解放军理工大学 | 基于Associated Hermite 正交函数的无条件稳定FDTD算法 |
CN104050147A (zh) * | 2013-03-13 | 2014-09-17 | 刘湘辉 | 将时域信号转换成频域信号的方法与系统 |
CN104375975A (zh) * | 2014-12-01 | 2015-02-25 | 天津工业大学 | 基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20100076732A1 (en) * | 2008-09-23 | 2010-03-25 | Sangpil Yoon | Meshfree Algorithm for Level Set Evolution |
WO2012109465A2 (en) * | 2011-02-09 | 2012-08-16 | Utilidata Inc. | Mesh delivery system |
-
2015
- 2015-11-06 CN CN201510751170.8A patent/CN105426339B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101944144A (zh) * | 2010-08-30 | 2011-01-12 | 陈玉君 | 一种基于无网格的布类仿真方法 |
CN102262699A (zh) * | 2011-07-27 | 2011-11-30 | 华北水利水电学院 | 基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法 |
CN104050147A (zh) * | 2013-03-13 | 2014-09-17 | 刘湘辉 | 将时域信号转换成频域信号的方法与系统 |
CN103970717A (zh) * | 2014-05-08 | 2014-08-06 | 中国人民解放军理工大学 | 基于Associated Hermite 正交函数的无条件稳定FDTD算法 |
CN104375975A (zh) * | 2014-12-01 | 2015-02-25 | 天津工业大学 | 基于双线性变换的一维真空Crank-Nicolson完全匹配层实现算法 |
Non-Patent Citations (3)
Title |
---|
地球物理正演中的无网格法算法研究;苏洲;《中国优秀硕士学位论文全文数据库 基础科学辑》;20130115(第01期);正文第二章和第三章 * |
基于无单元Galerkin法探地雷达正演模拟;冯德山 等;《地球物理学报》;20130131;第56卷(第1期);正文第4-5节 * |
基于混合边界条件的有限单元法GPR正演模拟;冯德山 等;《地球物理学报》;20121130;第55卷(第11期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN105426339A (zh) | 2016-03-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105426339B (zh) | 一种基于无网格法的线源时域电磁响应数值计算方法 | |
Dey et al. | Resistivity modelling for arbitrarily shaped two‐dimensional structures | |
Siripunvaraporn et al. | WSINV3DMT: vertical magnetic field transfer function inversion and parallel implementation | |
CN102798898B (zh) | 大地电磁场非线性共轭梯度三维反演方法 | |
CN105893678A (zh) | 一种时域有限差分的三维感应-极化双场数值模拟方法 | |
Yin et al. | 3D time-domain airborne EM forward modeling with topography | |
Blome et al. | Advances in three-dimensional geoelectric forward solver techniques | |
HU et al. | Pseudo‐three‐dimensional magnetotelluric inversion using nonlinear conjugate gradients | |
CN106980736A (zh) | 一种各向异性介质的海洋可控源电磁法有限元正演方法 | |
CN110852025B (zh) | 一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法 | |
Tselentis | An attempt to define Curie point depths in Greece from aeromagnetic and heat flow data | |
CN112949134A (zh) | 基于非结构有限元方法的地-井瞬变电磁反演方法 | |
Pan et al. | 2.5-D and 3-D DC resistivity modelling using an extrapolation cascadic multigrid method | |
CN108984939A (zh) | 基于3D Gauss-FFT的三维重力场正演方法 | |
CN112163332A (zh) | 基于自然单元无限元耦合的2.5维自然电场数值模拟方法 | |
Cai et al. | Effective 3D-transient electromagnetic inversion using finite-element method with a parallel direct solver | |
Persova et al. | Geometric 3-D inversion of airborne time-domain electromagnetic data with applications to kimberlite pipes prospecting in a complex medium | |
Xiao et al. | Fast 2.5 D and 3D inversion of transient electromagnetic surveys using the octree-based finite-element method | |
Zhang et al. | 3D inversion of large-scale frequency-domain airborne electromagnetic data using unstructured local mesh | |
Scheunert et al. | A cut-&-paste strategy for the 3-D inversion of helicopter-borne electromagnetic data—I. 3-D inversion using the explicit Jacobian and a tensor-based formulation | |
Liu et al. | Numerical modeling of the 2D time-domain transient electromagnetic secondary field of the line source of the current excitation | |
Wang et al. | Efficient 2D modeling of magnetic anomalies using NUFFT in the Fourier domain | |
CN114488327B (zh) | 基于地面基点的水平磁场与井中垂直磁场联合测量方法 | |
CN114114438B (zh) | 一种回线源地空瞬变电磁数据的拟三维反演方法 | |
CN114065577A (zh) | 一种直流电阻率小波伽辽金三维正演方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |