CN111521112B - 一种傅里叶及窗口傅里叶变换的联合相位重构算法 - Google Patents

一种傅里叶及窗口傅里叶变换的联合相位重构算法 Download PDF

Info

Publication number
CN111521112B
CN111521112B CN202010327645.1A CN202010327645A CN111521112B CN 111521112 B CN111521112 B CN 111521112B CN 202010327645 A CN202010327645 A CN 202010327645A CN 111521112 B CN111521112 B CN 111521112B
Authority
CN
China
Prior art keywords
window
phase
dimensional
distribution
fourier transform
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
CN202010327645.1A
Other languages
English (en)
Other versions
CN111521112A (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.)
Xian Technological University
Original Assignee
Xian Technological 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 Xian Technological University filed Critical Xian Technological University
Priority to CN202010327645.1A priority Critical patent/CN111521112B/zh
Publication of CN111521112A publication Critical patent/CN111521112A/zh
Application granted granted Critical
Publication of CN111521112B publication Critical patent/CN111521112B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B9/00—Measuring instruments characterised by the use of optical techniques
    • G01B9/02—Interferometers
    • G01B9/02083—Interferometers characterised by particular signal processing and presentation
    • G01B9/02084—Processing in the Fourier or frequency domain when not imaged in the frequency domain
    • 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/141—Discrete Fourier transforms
    • 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)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Discrete Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Signal Processing (AREA)
  • Instruments For Measurement Of Length By Optical Means (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

本发明公开了一种傅里叶及窗口傅里叶变换的联合相位重构算法,其特征在于,针对单帧载频干涉条纹图，利用傅里叶空频分析法滤除背景光强及其他无关项的干扰，利用傅里叶逆变换处理将其转换到空域内；随后对其进行二维窗口傅里叶变换处理，将得到的窗口傅里叶频谱分布中的最大值进行提取，滑动窗口函数获得整幅干涉条纹的窗口傅里叶脊分布，包裹相位可对窗口傅里叶脊求相位角计算求取。随后对该包裹相位进行二维解包裹处理，获得所需连续相位。本发明与传统WFT相比，克服了背景光强干扰问题，并通过正一级频谱移中处理化简了载频条纹在面形拟合时存在的消倾斜等繁杂过程，显著提高了WFT相位重构技术的精度，可进一步应用于大口径光学元件测量中。

Description

一种傅里叶及窗口傅里叶变换的联合相位重构算法
技术领域
本发明涉及光学干涉测量领域，尤其涉及一种傅里叶及窗口傅里叶变换的联合相位重构算法。
背景技术
光学干涉测量作为一种非接触测量方法被定义为新世纪精密测量技术的重要发展方向。相位提取是光学干涉测量中的关键一步。传统的相位提取技术主要分相移法及傅里叶变换法，但其各有局限性，相移法在测量大口径光学元件时易受到大气扰动，且处理数据较为复杂；傅里叶变换法有效的克服以上相移法存在的机械振动和气流影响，但傅里叶变换是一种全局分析技术，使局域内噪声滤除不彻底，降低相位重构精度。近些年来提出的窗口傅里叶变换相位提取算法在一定程度上解决了局部分析的问题，并通过局域(窗口)分析使干涉条纹不受噪声影响，成为当前大口径光学元件测量常用方法之一。
利用窗口傅里叶相位提取技术对某一干涉条纹进行处理时，默认窗口内包含的条纹成分为线性分布，但实际的载频条纹呈非线性分布，导致算法在运行时产生线性相位误差。目前对于干涉条纹的这种线性相位误差处理，一般是采取缩小窗口尺寸，使窗口内的干涉条纹近似认为线性，但小尺寸窗口在频域内尺寸激增，导致区域内条纹低频成分受到背景光强干扰，且小尺寸窗口处理下，算法运行时间冗长，不利于实际应用。
因此，提供一种新的傅里叶及窗口傅里叶变换的联合相位重构算法是本领域技术人员亟需解决的技术问题。
发明内容
为解决上述技术问题，本发明提供一种傅里叶及窗口傅里叶变换的联合相位重构算法
本发明提供一种傅里叶及窗口傅里叶变换的联合相位重构算法，包括以下步骤：
步骤1：对采集或仿真的初始载频条纹进行二维傅里叶变换，获得其频谱分布；
所述载频条纹强度表示为：
Figure BDA0002463796290000021
其中，
Figure BDA0002463796290000022
Id(x,y)是干涉条纹的背景光强分布，b(x,y)是干涉条纹的调制度分布，
Figure BDA0002463796290000023
为含有待测波面相位信息的相位分布函数，fx、fy为沿x,y方向上的空间载频，*表示复共轭；
步骤2：针对上述单帧载频条纹，对其进行二维傅里叶变换处理，由于载频成分的引入，使原始指数相位场内的正负一级频谱发生位置的偏移，在其对应频域内呈现正一,负一级频谱分离的状态,在此状态下，载频条纹对应频域分布可表示为：
F(f1,f2)＝A(f1,f2)+C(f1-fx,f2-fy)+C*(f1+fx,f2+fy) (2)
其中，等号右侧从左至右依次为：直流分量，正一级频谱及负一级频谱，载频的大小控制着正负一级频谱间的间距。
步骤3：根据不同条纹图的正一级频谱分布特点，选取合适的滤波窗类型，常用的几种滤波窗口，如汉宁窗，海明窗及布莱克曼滤波窗，可自行根据对应正一级频谱的分布特点进行选取；其次选定对应滤波窗口的中心带宽，其数值与正一级频谱的中心位置(fx,fy)保持一致，反复调试提取出完整的正一级频谱并滤除背景光强及无关项.将正一级频谱成分移至中心位置，获得滤除背景光强项的干涉条纹成分,记为C(f1,f2)。
步骤4：将步骤3提取的正一级频谱进行傅里叶逆变换处理，将其由频域转移到空域内，获得包含相位信息分布成分，记为c(x,y),表示为
Figure BDA0002463796290000031
步骤5：对步骤4获得的包含相位信息的正一级频谱的空域成分进行窗口傅里叶变换相位提取处理，对于WFT相位提取技术，输入端为上述c(x,y)，则其对应窗口傅里叶变换可表示为
Figure BDA0002463796290000032
根据仿真及实验的多次调试，选取WFT算法所用默认窗口类型为高斯窗，其一维表达式为
Figure BDA0002463796290000033
式中σ代表高斯窗函数在x轴方向上的延伸，可理解为高斯窗口尺寸大小.
对式(5)的一维窗口函数进行旋转或点乘，得到其对应二维分布，对应二维高斯窗函数可表示为
g(x,y)＝gx(x)gy(y) (6)
其中，gx(x)＝gy(y)＝g(x),定义σx、σy为此二维窗口函数在x,y轴方向上的延伸。在前期对窗口尺寸的精度分析中,确定了可使条纹线性相位误差抑制度达最大时的最优高斯窗口尺寸为σx＝σy＝5pixel；
将所述c(x,y)及二维高斯窗函数g(x,y)带入式(4)内，整理可得:Sf0(u,v；ξ,η)＝c(x,y)Gx,y(u,v；ξ,η) (7)
Gx,y(u,v；ξ,η)＝Gx,y(u,v；ξ)Gx,y(u,v；η) (8)
其中
Figure BDA0002463796290000041
Figure BDA0002463796290000042
定义Gx,y(u,v；ξ)、Gx,y(u,v；η)分别为沿x,y轴的一维复增益因子，而Gx,y(u,v；ξ,η)为二维复增益因子。
式(9)、(10)中的
Figure BDA0002463796290000043
分别表示在像素点(u,v)处沿x，y轴的局部频率，其在数值上等于相位的一阶偏导数；
Figure BDA0002463796290000044
分别表示在像素点(u,v)处沿x，y轴的局部曲率，其在数值上等于相位的二阶偏导数；
当在满足(ξ,η)＝(ωx,ωy)时，二维增益因子Gx,y(u,v；ξ,η)取得最大值，为
Figure BDA0002463796290000045
确定窗口覆盖区域内的最大谱值(或称为脊值)处的自变量值，将其赋值为该窗口内的局部频率(或称相位导数)，滑动窗口累计所有窗口内的脊值并叠加，获得整幅载频条纹的窗口傅里叶脊值分布，对其进行相位角计算获得包裹相位。
Figure BDA0002463796290000051
步骤:6：利用二维离散余弦相位解包裹算法对包裹相位进行解包处理获得连续相位，并利用Zernike多项式重构三维面形分布。
与相关技术相比，本发明具有以下有益效果：
1、传统WFT技术中所提取的包裹相位均包含载频成分，在后续面形拟合时还需人为进行消倾斜处理，本发明在对干涉条纹进行背景光强滤除后，将正一级频谱移中后消除了载频的影响，使后续包裹相位不包含载频，简化算法操作流程。
2、本发明解决了WFT相位提取技术在选取小尺寸窗口情况下出现的背景光强干扰问题，并通过最优窗口尺寸的共同作用对载频干涉条纹的线性相位误差做统一抑制，使该FT-WFT联合算法处理下的相位提取精度优于传统WFT技术。
附图说明
图1：基于本发明提出的FT-WFT联合相位重构算法流程图；
图2：经计算机仿真生成的单帧载频干涉条纹图；
图3：经傅里叶变换的载频干涉条纹频谱分布图；
图4：经矩形滤波窗提取的正一级频谱分布图
图5：移中的正一级频谱分布图；
图6：正一级频谱经傅里叶逆变换转换到空域内的强度分布图；
图7：滤除背景光强得到的干涉条纹成分经WFT后得到的窗口傅里叶脊值分布图
图8：对窗口傅里叶脊值求取相位角得到的包裹相位分布图；
图9：生成的包裹相位经二维离散余弦相位解包裹处理得到的连续相位分布图；
图10：传统WFT相位重构算法处理图2所获包裹相位分布图；
图11：包裹相位经Zernike多项式拟合重构的三维面形分布图；
具体实施方式
下面将结合附图和实施方式对本发明作进一步说明。
一种傅里叶及窗口傅里叶变换的联合相位重构算法，如图1所示，其详细步骤为：
步骤1：对单帧载频干涉条纹图进行二维傅里叶变换得到频谱分布；仿真尺寸为300×300pixel的载频干涉条纹图，添加均值μ＝0，方差σ＝1的高斯白噪声，获得初始干涉条纹图，如图2所示；对其进行二维傅里叶变换，得到其二维频谱分布，如图3所示为其对应三维频谱立体分布俯视图；
步骤2：正一级频谱的提取及转移；根据其正一级频谱分布特点，选取滤波窗口为矩形窗，并通过计算机调试，确定矩形窗口大小为24pixel，在频域内将干涉条纹傅里叶谱与滤波窗进行点乘运算，获得滤除背景光强及无关项的基频频谱，如图4所示，并将其移至区域中心位置，如图5所示；
步骤3：二维傅里叶逆变换得到空域分布；对移中的正一级频谱进行二维傅里叶逆变换处理，得到其空域内强度分布图，如图6所示；
步骤4：对IFT得到的二维条纹成分进行二维傅里叶变换得到窗口傅里叶谱；对经频谱处理后的载频条纹信号进行二维窗口傅里叶变换，其中窗口类型选取为高斯窗，窗口尺寸大小为5pixel，得到信号的窗口傅里叶谱，并通过最大谱值搜寻找到可使窗口区域内谱值达最大时的自变量取值，记为该点局部频率，图7所示为其脊值分布；对脊值进行相位角计算得到包裹相位，如图8所示；
步骤5：提取每个窗口内脊值(最大谱值)并对其做相位角计算获得包裹相位；利用离散余弦相位解包裹算法对步骤4获得的包裹相位进行解包处理，得到连续相位，如图9所示；
步骤6：进行离散余弦解包裹处理获得连续相位并重构面形；对初始干涉条纹图直接进行WFT技术处理，选用小尺寸窗口σ＝5pixel进行相位提取，所获包裹相位分布如图10所示；与本发明提取的包裹相位进行对比可知，利用传统WFT技术所获包裹相位存在相位在某些区域变化不明显呈连续现象，认为其部分低频成分受到背景光强干扰，本发明提出FT-WFT联合算法在对初始干涉条纹进行傅里叶空频处理，去除背景光强及其他无关项，所获包裹相位清晰，且不包含载频成分，无需做后续消倾斜处理；利用Zernike多项式重构面形得到图11所示的三维面形分布。
传统WFT技术中所提取的包裹相位均包含载频成分，在后续面形拟合时还需人为进行消倾斜处理，本发明在对干涉条纹进行背景光强滤除后，将正一级频谱移中后消除了载频的影响，使后续包裹相位不包含载频，简化算法操作流程。
本发明解决了WFT相位提取技术在选取小尺寸窗口情况下出现的背景光强干扰问题，并通过最优窗口尺寸的共同作用对载频干涉条纹的线性相位误差做统一抑制，使该FT-WFT联合算法处理下的相位提取精度优于传统WFT技术。
以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其它相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (1)

1.一种傅里叶及窗口傅里叶变换的联合相位重构算法,其特征在于：所述算法包括以下步骤：
步骤一：对采集或仿真的初始载频条纹进行二维傅里叶变换，获得其频谱分布；
所述载频条纹强度表示为：
Figure FDA0002976190500000011
其中，
Figure FDA0002976190500000012
Id(x,y)是干涉条纹的背景光强分布，b(x,y)是干涉条纹的调制度分布，
Figure FDA0002976190500000013
为含有待测波面相位信息的相位分布函数，fx、fy为沿x,y方向上的空间载频，*表示复共轭；
步骤二：针对单帧载频条纹，对其进行二维傅里叶变换处理，载频条纹对应频域分布表示为：
F(f1,f2)＝A(f1,f2)+C(f1-fx,f2-fy)+C*(f1+fx,f2+fy) (2)
其中，等号右侧从左至右依次为：直流分量，正一级频谱及负一级频谱，载频的大小控制着正负一级频谱间的间距；
步骤三：根据不同条纹图的正一级频谱分布特点，选取合适的滤波窗类型；其次选定对应滤波窗口的中心带宽，其数值与正一级频谱的中心位置(fx,fy)保持一致，反复调试提取出完整的正一级频谱并滤除背景光强及无关项,将正一级频谱成分移至中心位置，获得滤除背景光强项的干涉条纹成分,记为C(f1,f2)；
步骤四：将步骤三提取的正一级频谱进行傅里叶逆变换处理，将其由频域转移到空域内，获得包含相位信息分布成分，记为c1(x,y),表示为：
Figure FDA0002976190500000021
步骤五：对步骤四获得的包含相位信息的正一级频谱的空域成分进行窗口傅里叶变换相位提取处理，对于窗口傅里叶变换相位提取技术，输入端为上述c1(x,y)，其对应窗口傅里叶变换表示为
Figure FDA0002976190500000022
根据仿真及实验的多次调试，选取窗口傅里叶变换算法所用默认窗口类型为高斯窗，其一维表达式为
Figure FDA0002976190500000023
式中σ代表高斯窗函数在x轴方向上的延伸，可理解为高斯窗口尺寸大小；
对式(5)的一维窗口函数进行旋转或点乘，得到其对应二维分布，对应二维高斯窗函数可表示为
g(x,y)＝gx(x)gy(y) (6)
其中，gx(x)＝gy(y)＝g(x),定义σx、σy为此二维窗口函数在x,y轴方向上的延伸；在前期对窗口尺寸的精度分析中,确定了可使条纹线性相位误差抑制度达最大时的最优高斯窗口尺寸为σx＝σy＝5pixel；
将所述c1(x,y)及二维高斯窗函数g(x,y)带入式(4)内，整理可得
Sf0(u,v；ξ,η)＝c1(x,y)Gx,y(u,v；ξ,η) (7)
Gx,y(u,v；ξ,η)＝Gx,y(u,v；ξ)Gx,y(u,v；η) (8)
其中
Figure FDA0002976190500000031
Figure FDA0002976190500000032
定义Gx,y(u,v；ξ)、Gx,y(u,v；η)分别为沿x,y轴的一维复增益因子，而Gx,y(u,v；ξ,η)为二维复增益因子；
式(9)、(10)中的
Figure FDA0002976190500000033
分别表示在像素点(u,v)处沿x，y轴的局部频率，其在数值上等于相位的一阶偏导数；
Figure FDA0002976190500000034
分别表示在像素点(u,v)处沿x，y轴的局部曲率，其在数值上等于相位的二阶偏导数；
当在满足(ξ,η)＝(ωx,ωy)时，二维增益因子Gx,y(u,v；ξ,η)取得最大值，为
Figure FDA0002976190500000035
确定窗口覆盖区域内的最大谱值或称为脊值处的自变量值，将其赋值为该窗口内的局部频率或称相位导数，滑动窗口累计所有窗口内的脊值并叠加，获得整幅载频条纹的窗口傅里叶脊值分布，对其进行相位角计算获得包裹相位：
Figure FDA0002976190500000041
步骤六：利用二维离散余弦相位解包裹算法对包裹相位进行解包处理获得连续相位，并利用Zernike多项式重构三维面形分布。
CN202010327645.1A 2020-04-23 2020-04-23 一种傅里叶及窗口傅里叶变换的联合相位重构算法 Active CN111521112B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010327645.1A CN111521112B (zh) 2020-04-23 2020-04-23 一种傅里叶及窗口傅里叶变换的联合相位重构算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010327645.1A CN111521112B (zh) 2020-04-23 2020-04-23 一种傅里叶及窗口傅里叶变换的联合相位重构算法

Publications (2)

Publication Number Publication Date
CN111521112A CN111521112A (zh) 2020-08-11
CN111521112B true CN111521112B (zh) 2021-04-27

Family

ID=71904294

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010327645.1A Active CN111521112B (zh) 2020-04-23 2020-04-23 一种傅里叶及窗口傅里叶变换的联合相位重构算法

Country Status (1)

Country Link
CN (1) CN111521112B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114764149B (zh) * 2021-01-13 2023-04-07 中国石油化工股份有限公司 陡坡带砂砾岩体有利相带描述方法
CN115248549B (zh) * 2022-01-12 2024-05-24 浙江理工大学 一种打散消除杂散频谱噪声的数字全息三维重建方法
CN114964014B (zh) * 2022-05-23 2024-07-23 天津大学四川创新研究院 基于太赫兹的多层膜厚测量方法及可读存储介质
CN115307746B (zh) * 2022-08-10 2023-08-18 苏州科技大学 一种单色光的干涉测量相位计算方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5538936B2 (ja) * 2010-02-10 2014-07-02 キヤノン株式会社 解析方法、プログラム、記憶媒体、ｘ線位相イメージング装置
CN102628676B (zh) * 2012-01-19 2014-05-07 东南大学 一种光学三维测量中的自适应窗口傅里叶相位提取法
CN102620685B (zh) * 2012-03-23 2014-11-26 东南大学 一种基于史托克维尔变换的改进窗口傅里叶三维测量法
CN103267496B (zh) * 2013-05-20 2016-01-27 东南大学 一种基于小波变换的改进窗口傅里叶三维测量法
CN110160459B (zh) * 2019-06-06 2021-02-23 山东师范大学 基于改进光流场理论提取散斑干涉条纹图形变相位的方法

Also Published As

Publication number Publication date
CN111521112A (zh) 2020-08-11

Similar Documents

Publication Publication Date Title
CN111521112B (zh) 一种傅里叶及窗口傅里叶变换的联合相位重构算法
Kaufmann et al. Speckle noise reduction in television holography fringes using wavelet thresholding
Green et al. Investigation of the Fourier-transform method of fringe pattern analysis
CN107917676B (zh) 一种基于条纹图像频谱分析的干涉测量方法
Servin et al. Modern fringe pattern analysis in interferometry
CN102628676A (zh) 一种光学三维测量中的自适应窗口傅里叶相位提取法
CN103267496B (zh) 一种基于小波变换的改进窗口傅里叶三维测量法
CN104155011B (zh) 一种二维干涉图的相位提取方法
CN107167087A (zh) 一种基于经验小波变换的物体变形测量方法
CN106931905B (zh) 一种基于非线性优化的数字莫尔条纹相位提取方法
CN112381731B (zh) 一种基于图像去噪的单帧条纹图像相位分析方法及***
Gdeisat et al. Fast fringe pattern phase demodulation using FIR Hilbert transformers
CN104614083B (zh) 一种恢复相移干涉图相位分布、及获取两幅图间相移量的方法
Trusiak et al. Space domain interpretation of incoherent moiré superimpositions using FABEMD
Li et al. Quadratic polar coordinate transform technique for the demodulation of circular carrier interferogram
CN116224327B (zh) 一种基于学习网络的矿区大梯度形变区相位解缠方法
CN105404129B (zh) 基于内积算法的三步任意相移消除数字全息零级像的方法
Cywińska et al. Enhancing single-shot fringe pattern phase demodulation using advanced variational image decomposition
CN110500968B (zh) 基于稀疏傅里叶变换的数字莫尔干涉相位实时测量方法
Kulkarni et al. Patch-wise denoising of phase fringe patterns based on matrix enhancement
Abdul-Rahman et al. Three-dimensional Fourier fringe analysis
Kujawinska et al. Fringe pattern analysis using Fourier transform techniques
CN115307746B (zh) 一种单色光的干涉测量相位计算方法
Montrésor et al. Processing of digital holograms with 2d windowed Fourier transform: Study, performance analysis and gpu implementation
CN114623782B (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