CN109633527B - 基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法 - Google Patents
基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法 Download PDFInfo
- Publication number
- CN109633527B CN109633527B CN201811531613.2A CN201811531613A CN109633527B CN 109633527 B CN109633527 B CN 109633527B CN 201811531613 A CN201811531613 A CN 201811531613A CN 109633527 B CN109633527 B CN 109633527B
- Authority
- CN
- China
- Prior art keywords
- time difference
- microphone
- arrival time
- arrival
- matrix
- 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
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
- G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
- G01S3/802—Systems for determining direction or deviation from predetermined direction
- G01S3/8027—By vectorial composition of signals received by plural, differently-oriented transducers
-
- Y—GENERAL 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
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Abstract
本发明公开了一种基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,该方法包括以下步骤:首先根据嵌入式平面麦克风阵列获取的实测数据进行波达时间差估计;之后利用波达时间差估计值,根据波达时间的大小关系选择通视麦克风作为参考麦克风;然后基于选取的参考麦克风构建波达时间差矩阵,进行施加低秩及几何约束条件下的凸优化校正;最后根据校正后的波达时间差结果构建波达时间差向量,进而进行最小二乘方位角估计。本发明显著减小了嵌入式平面麦克风阵列最小二乘方位角估计偏差的统计均值与方差,具有广阔的应用前景。
Description
技术领域
本发明属于嵌入式麦克风阵列声源测向领域,特别是一种基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法。
背景技术
麦克风阵列声源测向方法很多,包括高分辨谱估计、波束形成、波达时间差(TDOA)等。其中,麦克风间TDOA测量结合最小二乘(LS)估计的两级声源测向方法因实现复杂度低、环境适应性好等优点在各种实时声探测任务中得到广泛应用。因受外形设计、空间尺寸、使用运输等多种条件限制,图3所示的便携式枪声探测系统所用的麦克风阵列通常采用嵌入式安装方式,可能导致由于外壳遮挡而造成目标与某个或某些麦克风之间处于非通视(NLOS)状态。考虑到经过远距离传播后的枪声膛口波信号能量主要集中在数百Hz的低频范围,非通视麦克风受壳体衍射传播效应影响,相关麦克风间的实际TDOA不再与麦克风间距和目标位置保持理想全通视条件下的理论对应关系,进而造成基于LS测向方法的膛口波方向估计结果出现明显偏差。
有关非通视情况下的声源目标测向问题,许多文献主要关注声探测系统整体因障碍阻隔而无法直视声源目标时的解决方案,一般不涉及壳体遮挡对嵌入式阵列目标方向估计的影响及其处理方法。针对仅有部分非通视麦克风的情形,已有研究表明,采用球麦克风阵列声场模型能有效抑制衍射传播干扰、获得精确的声源目标方位估计。但基于模型的衍射传播效应抑制方法一般较为复杂,难以实时实现;同时受成本、结构、使用方式等制约,便携式枪声探测系统通常采用麦克风数很少的平面阵结构,很难严格满足均匀球形阵列甚至圆形阵列的模型要求。
发明内容
本发明所要解决的技术问题在于针对嵌入式平面麦克风阵列的声源测向误差,提供一种实时性高、易于实现、且测向较为精确的嵌入式平面麦克风阵列声源测向方法。
实现本发明目的的技术解决方案为:基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,包括以下步骤:
步骤1、通过嵌入式平面麦克风阵列获取实测数据,并对实测数据进行波达时间差估计;
步骤2、利用步骤1中波达时间差估计值,根据波达时间的大小关系选择通视麦克风作为参考麦克风;
步骤3、基于步骤2中选择的参考麦克风构建波达时间差矩阵,进行施加低秩及几何约束条件下的凸优化校正;
步骤4、根据步骤3中校正后的波达时间差结果构建波达时间差向量,进而进行最小二乘方位角估计,完成嵌入式平面麦克风阵列声源测向。
本发明与现有技术相比,其显著优点为:1)本发明中通过选择距离声源最近的麦克风、施加低秩以及几何约束,降低了测向误差的均值及标准差,提高了测向精度及稳健性;2)本发明中选择距离声源最近的麦克风,尽可能避免了壳体遮挡对测向性能的影响;3)本发明采用反对称TDOA矩阵的低秩约束条件,抑制了实测TDOA的测量随机误差;4)本发明采用阵列形状的几何约束条件,抑制了壳体遮挡引起的衍射传播误差;5)本发明中将TDOA的校正转化为带约束条件的凸优化问题,利用拉格朗日乘子法求解凸优化问题的解析解,使得该方法的计算量较小,可满足系统的实时性需要,也降低了成本;6)本发明中通过将TDOA校正转化为带约束条件的凸优化问题,仅通过改变各矩阵的维数和几何约束的表现形式,就可针对性的推广到麦克风数量或阵列形状不同的其它麦克风阵列系统。
下面结合附图对本发明作进一步详细的描述。
附图说明
图1为本发明基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法流程图。
图2为本发明实施例中实验布设示意图。
图3为本发明实施例所用的嵌入式平面麦克风阵列示意图。
图4为本发明实施例中参考麦克风选择的结构框图。
图5为本发明方法和常规方法的方位角估计误差分布对比示意图。
图6为本发明方法和常规方法的方位角估计误差分布直方图,其中(a)为常规方法的方位角估计误差分布直方图,(b)为本发明方法的方位角估计误差分布直方图。
图7为本实施例中均值补偿后两种方法的绝对方位角估计误差累积概率分布示意图。
具体实施方式
结合图1,本发明基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,包括以下步骤:
步骤1、通过嵌入式平面麦克风阵列获取实测数据,并对实测数据进行波达时间差估计。
步骤2、利用步骤1中波达时间差估计值,根据波达时间的大小关系选择通视麦克风作为参考麦克风。
进一步地,步骤2具体为:
步骤2-1、根据波达时间差估计值构建波达时间差估计矩阵,具体为:
令远场目标信号传播至麦克风Mi的波达时间为ti,i=1~n;定义嵌入式平面麦克风阵列的波达时间向量为x=[t1,t2,t3,…,tn]T,远场目标信号传播至麦克风Mi和Mj之间的波达时间差τij=ti-tj,结合维度为n的全1向量1=[1,1,1,…,1]T可得相应的反对称波达时间差矩阵M:
M=x·1T-1·xT
步骤2-3、以记录最多的行序号即众数作为距离目标最近的通视麦克风的序号,并选取该麦克风作为参考麦克风。
进一步地,n的取值为4。
步骤3、基于步骤2中选择的参考麦克风构建波达时间差矩阵,进行施加低秩及几何约束条件下的凸优化校正。
进一步地,步骤3具体为:
施加波达时间差矩阵代数特征的低秩约束,以及对应阵列形状特征的几何约束,将对麦克风全通视情况下的理想波达时间差矩阵M的估计任务转化为有约束凸优化问题:
式中,||·||F为F-范数,约束矩阵D=[1,g],包括低秩约束1Tx=0和几何约束gTx=0,前者与理想全通视TDOA矩阵M应为秩2的低秩要求相关,用于抑制高斯随机测量噪声,后者则体现了平行四边形阵列所对应的形状特征几何约束,其中g称为约束向量;
对上式应用拉格朗日乘数法构造函数:
式中,待定系数λ1、λ2为拉格朗日乘子;
考虑:
分别对x、λ1、λ2的一阶导数为0,可得
步骤4、根据步骤3中校正后的波达时间差结果构建波达时间差向量,从而完成最小二乘方位角估计。
进一步地,步骤4具体为:
步骤4-1、构建目标方向对应的麦克风间波达时间差向量b为:
式中,A为与嵌入式平面麦克风阵列尺寸及形状有关的阵列系数矩阵,为远场目标方向矢量,θ、分别为远场目标方位角和俯仰角,τi1为第i个麦克风的坐标向量与第1个麦克风的坐标向量的差值,mi为第i个麦克风的坐标向量,c为声速;
步骤4-2、求取上式向量b的线性最小二乘解为:
步骤4-3、进行最小二乘方位角估计为:
下面结合实施例对本发明作进一步详细的描述。
实施例
结合图1,本实施例基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,包括以下内容:
步骤1、本实施例的实验在较空旷的室外环境进行,实验布设示意图如图2所示。实施例中采用的嵌入式平面麦克风阵列示意图如图3所示,其包括4个麦克风,嵌入式平面麦克风阵列系统被水平固定在三脚架上,以阵列中心为坐标原点,通过水平顺时针旋转阵列系统等效从小到大改变目标方位角。实施例中作为模拟目标的扬声器与嵌入式平面麦克风阵列高度基本相同,放置在距离阵列中心15m处。本实施例中,目标方位角从0°开始以30°步进递增变化一周,在每个目标方位角位置进行26次独立实验,12个方位角位置共计进行312次独立实验。每次实验中扬声器播放一段含有外场实测的枪声膛口波信号录音。对312次独立实验的实测数据各自进行波达时间差估计。
步骤2、对于一组波达时间差估计值,显然通视麦克风离声源最近,因此根据波达时间的大小关系选择通视麦克风,并作为参考麦克风。结合图4,具体如下:
步骤2-1、令嵌入式平面麦克风阵列系统中的其中一个麦克风为M1,其余麦克风以俯视逆时针为序依次命名M2、M3、M4,远场目标信号传播至麦克风Mi的波达时间为ti,i=1~4;定义阵列的波达时间向量为x=[t1,t2,t3,t4]T,远场目标信号传播至麦克风Mi和Mj之间的波达时间差τij=ti-tj,i、j=1~4,结合全1向量1=[1,1,1,1]T可得相应的反对称波达时间差矩阵M:
步骤2-3、以记录最多的行序号即Tmin中的众数作为距离目标最近的通视麦克风序号,并选取该麦克风作为参考麦克风。
步骤3、基于步骤2中选择的参考麦克风构建波达时间差矩阵,进行施加低秩及几何约束条件下的凸优化校正,具体如下:
施加波达时间差矩阵代数特征的低秩约束,以及对应阵列形状特征的几何约束,将对麦克风全通视情况下的理想波达时间差矩阵M的估计任务转化为有约束凸优化问题:
式中,||·||F为F-范数,约束矩阵D=[1,g],包括低秩约束1Tx=0和几何约束gTx=0,前者与理想全通视波达时间差矩阵M应为秩2的低秩要求相关,用于抑制高斯随机测量噪声,后者则体现了平行四边形阵列所对应的形状特征几何约束,由于τ12=τ43,因而式中g=[1,-1,1,-1]T。
针对上式应用Lagrange乘数法构造函数:
式中,待定系数λ1、λ2为拉格朗日乘子;
考虑:
分别对x、λ1、λ2的一阶导数为0,可得
步骤4、根据步骤3中校正后的波达时间差结果构建波达时间差向量,从而完成最小二乘方位角估计,具体如下:
步骤4-1、对于嵌入式阵列n=4,构建目标方向对应的麦克风间波达时间差向量b为:
式中,A为与阵列尺寸及形状有关的阵列系数矩阵,为远场目标方向矢量,θ和分别为远场目标方位角和俯仰角,τi1为第i个麦克风的坐标向量与第1个麦克风的坐标向量的差值,mi为第i个麦克风的坐标向量,c为声速;
步骤4-2、求取上式向量b的线性最小二乘解为:
步骤4-3、进行最小二乘方位角估计为:
步骤5、采用理想全通视直线传播模型下基于TDOA的常规LS测向算法(下文简称常规方法)作为性能对比,并定义方位角估计误差为目标方位角估计值与真实值之差。
图5为本发明方法和常规方法在12个测试方位角上各26次独立实验的方位角估计误差分布。图中显示,常规方法的方位角估计误差在一些方位角处散布较大,而经过本发明方法校正后,各方位角的方位角估计误差分布则已趋向集中。
图6比较了两种方法的方位角估计误差总体分布直方图。对比图6(a)、(b)可以看出,常规方法的方位角估计误差主要集中分布在–2.03°~5.14°范围内(上下四分位数误差范围),分布均匀程度较高,且误差绝对值的最大值高达23°。而本发明方法方位角估计误差则主要集中分布在–1.38°~2.84°范围内(上下四分位数误差范围),相应频度随着误差绝对值增大而很快下降,误差绝对值的最大值也降为10°。
此外,常规方法和本发明方法的方位角估计误差统计均值±方差分别为1.73±5.87°和0.90±3.55°,统计中值分别为1.88°和0.87°。实际目标方位角通常未知且随机,通常假设实际目标方位角在0°~360°之间均匀分布。为使全方位的总体测角性能相对均衡,同时考虑到两种方法方位角估计误差的统计均值与统计中值基本相等,实际应用时可预先补偿上述统计均值。两种方法在均值补偿后的方位角估计误差绝对值所对应的累积概率分布对比如图7所示。需要强调的是,就枪声探测系统实战环境下的技术指标而言,当前以英国神听猛击(Ears SWATS)II代、美国飞镖战警-X(Boomerang Warrior-X)系统为代表的国外同类便携式产品的枪声膛口波测向性能通常以±7.5°为测向精度考核标准。通过考察图7可知,针对方位角估计误差绝对值在7.5°以内的累积概率,常规方法仅为0.86,而与之形成鲜明对比的是,本发明方法达到0.94,测向性能提升明显。
综上所述,本发明的方法通过在低秩约束的基础上增加阵列形状几何约束来构建凸优化问题并选取距离目标最近的通视麦克风作为参考麦克风,能够同时有效抑制TDOA随机测量噪声与壳体遮挡衍射传播效应对基于TDOA的LS测向方法的不利影响,方位角估计精度得到明显提高。需要说明的是,仅通过改变各矩阵或向量的维数以及几何约束的表现形式,上述测向方法就可直接推广到麦克风数量或阵列形状各不相同的其它嵌入式麦克风阵列系统。
Claims (5)
1.基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,其特征在于,包括以下步骤:
步骤1、通过嵌入式平面麦克风阵列获取实测数据,并对实测数据进行波达时间差估计;
步骤2、利用步骤1中波达时间差估计值,根据波达时间的大小关系选择通视麦克风作为参考麦克风;
步骤3、基于步骤2中选择的参考麦克风构建波达时间差矩阵,进行施加低秩及几何约束条件下的凸优化校正;
步骤4、根据步骤3中校正后的波达时间差结果构建波达时间差向量,进而进行最小二乘方位角估计,完成嵌入式平面麦克风阵列声源测向。
2.根据权利要求1所述的基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,其特征在于,步骤2中所述利用步骤1中波达时间差估计值,根据波达时间的大小关系选择通视麦克风作为参考,具体如下:
步骤2-1、根据波达时间差估计值构建波达时间差估计矩阵,具体为:
令远场目标信号传播至麦克风Mi的波达时间为ti,i=1~n;定义嵌入式平面麦克风阵列的波达时间向量为x=[t1,t2,t3,…,tn]T,远场目标信号传播至麦克风Mi和Mj之间的波达时间差τij=ti-tj,结合维度为n的全1向量1=[1,1,1,…,1]T可得相应的反对称波达时间差矩阵M:
M=x·1T-1·xT
步骤2-3、以记录最多的行序号即众数作为距离目标最近的通视麦克风的序号,并选取该麦克风作为参考麦克风。
3.根据权利要求2所述的基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,其特征在于,步骤2-1中n=4。
4.根据权利要求3所述的基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,其特征在于,步骤3中所述基于步骤2中选择的参考麦克风构建波达时间差矩阵,进行施加低秩及几何约束条件下的凸优化校正,具体如下:
施加波达时间差矩阵代数特征的低秩约束,以及对应阵列形状特征的几何约束,将对麦克风全通视情况下的理想波达时间差矩阵M的估计任务转化为有约束凸优化问题:
式中,||·||F为F-范数,约束矩阵D=[1,g],包括低秩约束1Tx=0和几何约束gTx=0,其中g称为约束向量;
5.根据权利要求1所述的基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法,其特征在于,步骤4所述根据步骤3中校正后的波达时间差结果构建波达时间差向量,进而进行最小二乘方位角估计,具体为:
步骤4-1、构建目标方向对应的麦克风间波达时间差向量b为:
式中,A为与嵌入式平面麦克风阵列尺寸及形状有关的阵列系数矩阵,为远场目标方向矢量,θ、分别为远场目标方位角和俯仰角,τi1为第i个麦克风的坐标向量与第1个麦克风的坐标向量的差值,mi为第i个麦克风的坐标向量,c为声速;
步骤4-2、求取上式向量b的线性最小二乘解为:
步骤4-3、进行最小二乘方位角估计为:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811531613.2A CN109633527B (zh) | 2018-12-14 | 2018-12-14 | 基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811531613.2A CN109633527B (zh) | 2018-12-14 | 2018-12-14 | 基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109633527A CN109633527A (zh) | 2019-04-16 |
CN109633527B true CN109633527B (zh) | 2023-04-21 |
Family
ID=66073935
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811531613.2A Active CN109633527B (zh) | 2018-12-14 | 2018-12-14 | 基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109633527B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US11902755B2 (en) | 2019-11-12 | 2024-02-13 | Alibaba Group Holding Limited | Linear differential directional microphone array |
CN112714383B (zh) * | 2020-12-30 | 2022-03-11 | 西安讯飞超脑信息科技有限公司 | 麦克风阵列的设置方法、信号处理装置、系统及存储介质 |
CN116017281B (zh) * | 2022-12-30 | 2023-10-24 | 深圳市中承科技有限公司 | 一种基于超宽带通信技术的室内定位方法 |
Family Cites Families (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP3863323B2 (ja) * | 1999-08-03 | 2006-12-27 | 富士通株式会社 | マイクロホンアレイ装置 |
CN1235429C (zh) * | 2002-11-07 | 2006-01-04 | 华为技术有限公司 | 一种位置估计方法 |
JP5278365B2 (ja) * | 2010-03-31 | 2013-09-04 | 沖電気工業株式会社 | 位置推定装置及び位置推定方法 |
CN105068042B (zh) * | 2015-07-29 | 2017-10-13 | 中国科学院上海微系统与信息技术研究所 | 一种基于麦克风阵列的运动目标计数方法 |
US9989633B1 (en) * | 2017-03-15 | 2018-06-05 | Cypress Semiconductor Corporation | Estimating angle measurements for source tracking using a phased array system |
CN107132505B (zh) * | 2017-05-19 | 2019-11-08 | 中国人民解放军信息工程大学 | 直达与非直达混合场景中的多目标直接定位方法 |
CN107037397A (zh) * | 2017-06-21 | 2017-08-11 | 哈尔滨工业大学 | 一种在波达方向估计中校正多种阵列误差的方法 |
CN107770859A (zh) * | 2017-09-21 | 2018-03-06 | 天津大学 | 一种考虑基站位置误差的tdoa‑aoa定位方法 |
-
2018
- 2018-12-14 CN CN201811531613.2A patent/CN109633527B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN109633527A (zh) | 2019-04-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109633527B (zh) | 基于低秩及几何约束的嵌入式平面麦克风阵列声源测向方法 | |
US9949033B2 (en) | Planar sensor array | |
Dai et al. | A sparse representation method for DOA estimation with unknown mutual coupling | |
CN108549059B (zh) | 一种复杂地形条件下的低空目标仰角估计方法 | |
CN107919535B (zh) | 一种基于定向双圆阵的立体阵列天线及其构建方法 | |
CN107884741A (zh) | 一种多球阵列多宽带声源快速定向方法 | |
US7839721B1 (en) | Modal beam processing of acoustic vector sensor data | |
CN109581388B (zh) | 一种实时三维成像声纳的近场宽视角波束形成方法 | |
CN114325584B (zh) | 基于合成孔径的多阵元超声波声源三维成像方法及系统 | |
CN105445718A (zh) | 一种基于阵列重构的分布式多载舰超视距雷达的doa估计方法 | |
CN111025273A (zh) | 一种畸变拖曳阵线谱特征增强方法及系统 | |
CN107861096A (zh) | 基于声音信号到达时间差的最小二乘测向方法 | |
CN110736976B (zh) | 一种任意阵形的声纳波束形成器性能估计方法 | |
CN113835063B (zh) | 一种无人机阵列幅相误差与信号doa联合估计方法 | |
CN107202975B (zh) | 一种二维矢量阵阵元姿态误差校正方法 | |
US9578433B2 (en) | Method for self-calibration of a set of sensors, in particular microphones, and corresponding system | |
CN111983599B (zh) | 一种基于方位-俯仰字典的目标二维doa估计方法 | |
CN109541526A (zh) | 一种利用矩阵变换的圆环阵方位估计方法 | |
CN110554358B (zh) | 一种基于虚拟球阵列扩展技术的噪声源定位识别方法 | |
CN109669172B (zh) | 基于主瓣内强干扰抑制的弱目标方位估计方法 | |
CN111323746A (zh) | 一种双圆阵的方位-等效时延差被动定位方法 | |
CN114167346B (zh) | 基于协方差矩阵拟合阵元扩展的doa估计方法及系统 | |
CN111431575B (zh) | 基于常规波束形成的来波方向稀疏重构方法 | |
Jo et al. | Sine-based EB-ESPRIT for source localization | |
CN113075633A (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 |