CN102175989B - 一种非相干分布式信号二维波达角的测定方法 - Google Patents
一种非相干分布式信号二维波达角的测定方法 Download PDFInfo
- Publication number
- CN102175989B CN102175989B CN 201010618502 CN201010618502A CN102175989B CN 102175989 B CN102175989 B CN 102175989B CN 201010618502 CN201010618502 CN 201010618502 CN 201010618502 A CN201010618502 A CN 201010618502A CN 102175989 B CN102175989 B CN 102175989B
- Authority
- CN
- China
- Prior art keywords
- angle
- array
- battle array
- pitch
- circle battle
- 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.)
- Expired - Fee Related
Links
Images
Landscapes
- Mobile Radio Communication Systems (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
该发明属于无线移动通信技术领域中采用阵列天线测定非相干多径分布式信号二维波达角的测定方法。包括:设置由两个均匀圆阵组成的立体天线阵列,建立两个圆阵数据的协方差矩阵,确定子矩阵及中心俯仰角初值,建立任意一个圆阵的协方差矩阵,建立伪噪声子矩阵,确定二维谱函数,最后分别确定中心方位角及确定中心俯仰角。该发明对中心方位角和中心俯仰角测定的均方根误差仅为角度扩展的十分之一左右,并避免了多维谱搜索和非线性优化。因而,本发明具有方法简单、可靠,可大幅度降低计算处理量及处理的复杂程度,以及测定的效率和精确度高,应用于智能天线系统中、可提高智能天线系统在多径环境中的性能及实际应用的价值等特点。
Description
技术领域
本发明属于无线移动通信技术领域,特别是涉及一种采用阵列天线测定非相干多径分布式信号(简称非相干分布信号)二维波达方向的二维波达角(Direction of Arrival,简称DOA)的测定方法。
背景技术
随着近年来移动通信事业的飞速发展,通信业务需求量越来越大,特别是三代移动通信等概念的提出,对通信技术提出了更高的要求。在现代复杂的移动通信环境和频带资源限制下要达到极大的通信容量和极好的通信质量等目标,主要受到多个因素限制:由于移动环境的复杂性,通常从发射机到接收机的信号包含有折射、衍射、绕射等多种信号成分,由于多径效应的影响产生信号的时延扩展,不同多径分量的到达时间不同,造成接收信号的时域展宽,直接导致接收信号的码间干扰;在码分多址(CDMA)系统中,由于不同用户扩频码的非正交性或异步将造成信号的共信道干扰(CCI)。对于这些干扰,仅仅采用目前的数字通信技术是远远无法克服的。近来展开研究的移动通信中智能天线技术,为克服和减轻这些限制,提供了有力的技术支持。在无线移动通信系统中采用智能天线技术,实现最优收发天线波束形成和自适应信号处理,可为蜂窝系统提供高质量的数据链路,提高基站天线的覆盖范围以及系统容量和业务质量,降低移动用户的码间干扰和多址干扰,以及降低所需发射功率等。
在智能天线系统中,准确估计来波的到达方向(DOA)非常重要,它不但有助于智能天线上、下行链路中波束更精确地指向某个用户,更好地进行信息的定向发送和接收;而且可为空时信号处理中其它处理方法提供重要的技术支持。
无线通信系统的电波传播有直达波传播和非直达波传播两种方式。对于直达波传播,电波在接收机和发射机之间以直线波的方式传播,其条件是第一涅费尔区内无障碍物。在陆地环境下,起伏的地形、建筑物和植被等传播环境的影响,直达波传播非常困难。对于非直达波传播,接收和发射机之间的传播路径上存在障碍物而无直射路线,到达接收机的信号都是通过反射、衍射和散射等途径。相对于直达波而言,这几种传播途径均会造成信号的时延、衰耗,极化方式和稳定性都有很大的不同,严重影响传输质量。工程上,对于传播路线上的遮蔽物,可以采用增加发射功率、或者提高收发天线有效高度加以避开等措施。但是,功率的加大和高度的提升是有限度的,并且在CDMA系统中,往往会增加对周围小区的同频干扰,对系统容量极为小利。
在直达波传播情况下,移动台和观测阵列之间存在直接路线,因此一般可将信号源假设为点源。在非直达波传播条件下,当多径数目较大或者难以有效地分离各条路径时,传统的点源模型及其参数估计方法不再适用,此时,采用空间分布式信号模型来描述信号的多径传播特性更适合。由于采用分布式信号模型,只需要少数几个分布参数就能完整地描述非直达波的分布特性,而不用区分所有的传播路径,数学上处理起来非常方便。
分布源根据其源内不同波达方向之间是否相干,可分为相干分布(CoherentlyDistributed,CD)源和非相干分布(Incoherently Distributed,ID)源。当信道缓慢变化,或者信号遇到光滑的物体反射时,通常会产生相干分布信号;而当信道快速变化,或者当信号遇到粗糙的物体反射时,则会产生非相干分布信号。在移动通信中,非相干分布信号是普遍存在的;因此,研究非相干分布信号方位的测定方法具有十分重要的实际意义。
对于非相干分布信号的波达角测定问题,在过去的十多年中,人们已做了大量深入的研究,并提出了不少方法。例如,DSPE(distributed signal parameter estimator)、DISPARE(dispersed signal parametric estimation)方法,ESPRIT(rotational variance technique)方法,基于协方差匹配的COMET(covariance fitting approach)和COMET-EXIP方法,以及基于波束形成的广义SMVDR(spread minimum variance distortionless response)和广义Capon方法等。然而,上述方法大多为研究非相干分布源一维波达方向的测定问题,其中仅文献“BOUJEMAA H.Extension of COMET algorithm to multiple diffuse source localizationin azimuth and elevation[J].Euro Trans Telecomms,2005,16:557-566.”将COMET方法扩展应用于非相干分布源的二维波达方向的测定。该方法首先确定样本协方差矩阵,然后基于理想协方差矩阵和样本协方差矩阵之间的最小均方原则建立代价函数,最终由四维非线性最优化确定所测参数。COMET方法虽然检测精度高、接近于最大似然方法,但它涉及中心方位角、中心俯仰角、方位角扩展和俯仰角扩展的四维非线性优化问题,因而存在计算量巨大、处理复杂度很高,而实时性差,极不利于工程应用等缺陷。
发明内容
本发明的目的是设计一种非相干分布式信号二维波达角的测定方法,以达到降低数据处理的复杂度、简化处理程序,有效提高对非相干分布源二维波达角(DOA)测定的效率、精确度及实时性,应用于智能天线系统中、可提高智能天线系统在多径环境中的性能及实际应用的价值等目的。以克服背景技术或处理方法复杂度过高、或者不能应用于二维以及实时性差、难以实际应用等缺陷。
本发明的解决方案是针对背景技术存在的计算处理量巨大、复杂度高等缺陷,首先利用双平行圆阵获得各分布信号中心俯仰角的初值,再由阵列广义方向矢量和伪噪声子空间之间的正交性构造一种二维谱函数,然后通过二维谱函数并利用所得中心俯仰角的初值得到各分布信号的中心方位角的一维谱函数,最后通过一维搜索即可测出分布信号的中心方位角;再由测得的中心方位角与上述二维谱函数得到关于中心俯仰角的一维函数,通过一维搜索又进一步测出分布信号的中心俯仰角;从而大幅度降低处理量及处理的复杂程度、有效提高测定的效率和精确度,从而实现其发明目的。因此,本发明方法包括;
步骤1.设置立体天线阵列:在三维空间中设立一个由两个均匀圆阵(L1和L2)组成的立体天线阵列;每个圆阵的半径及阵元数(K)均相同,其中一个圆阵(L1或L2)设于x-y平面内、另一个圆阵则平行于该圆阵设置,以接收各非相干分布信号的数据矢量;
步骤2.建立两个圆阵数据的协方差矩阵:工作时首先将两圆阵(L1和L2)中各阵元每次接收的数据分别合并为一列、以建立协方差矩阵
步骤3.确定子矩阵:对步骤2所得协方差矩阵进行特征分解,并利用最大特征值所对应的特征向量建立特征向量矩阵(ES)、最大特征的个数为分布信号数的3倍,然后再将所得特征向量矩阵(ES)以阵元数为行数、将其分为前后两个子矩阵(ES1和ES2);
步骤4.确定中心俯仰角初值:利用总体最小二乘法(Total Least-Squares,TLS)建立步骤3所得两个子矩阵(ES1和ES2)的旋转矩阵(Ψ),并利用该旋转矩阵的特征值确定各分布信号中心俯仰角的初值;
步骤5.建立任意一个圆阵(L1或L2)的协方差矩阵:根据所选圆阵各次接收的数据建立该圆阵的协方差矩阵;
步骤6.确定伪噪声子矩阵:对步骤5得到的协方差矩阵进行特征分解、得特征值和特征向量,并根据单个圆阵的阵元数与分布信号个数的1.5-2.5倍值作为伪噪声子矩阵的维数,并由该维数个小特征值所对应的特征向量建立伪噪声子矩阵;
步骤7.确定二维谱函数:利用任意一个圆阵(L1或L2)的方向矢量,及该方向矢量对方位角的偏导数和对俯仰角的偏导数、建立广义方向矢量,并利用该广义方向矢量与伪噪声子矩阵之间的正交性确定二维谱函数;
步骤8.确定中心方位角:将步骤4所得各中心俯仰角的初值依次代入步骤7所得二维谱函数中,并用一维搜索方法确定出各初值所对应的中心方位角;
步骤9.确定中心俯仰角:将步骤8所得各中心方位角依次代入步骤7所得二维谱函数中,仍通过一维搜索方法确定出各中心方位角所对应的中心俯仰角。
所述每个圆阵的阵元数(K)均相同,其阵元数K=3,4,…20。
所述建立两个圆阵数据的协方差矩阵为:
其中:N为采样次数,n为采样的序号,z(n)H为z(n)的转置共轭,z(n)为由圆阵L1和L2的各次接收数据合并的一列,即:
在步骤3中所述特征向量矩阵为ES={e1,e2,…,e3q},式中,e1,e2,…,e3q为最大的(3q)个特征值所对应的特征向量;而所述两个子矩阵ES1和ES2则由ES按下式分为(前后)两个子矩阵
其中,q为分布信号的个数,K为单个圆阵的阵元数。
所述各分布信号中心俯仰角的初值由下式确定:
其中:i为各特征值组的序号、i=1,2,…q,3(i-1)+m为旋转矩阵Ψ特征值ξ的下标,m为各对应特征值组中各元素的序号。
所述任意一个圆阵(L1或L2)的协方差矩阵为:
其中:N为采样次数,n为采样的序号,x(n)H为x(n)的转置共轭。
所述二维谱函数为:
本发明由于首先利用双平行圆阵获得各分布信号中心俯仰角的初值,再由阵列广义方向矢量和伪噪声子空间之间的正交性构造一种二维谱函数,然后通过二维谱函数并利用所得中心俯仰角的初值得到各分布信号的中心方位角的一维谱函数,最后通过一维搜索即可测出分布信号的中心方位角;再由测得的中心方位角与上述二维谱函数得到关于中心俯仰角的一维函数,通过一维搜索又进一步测出分布信号的中心俯仰角;从而大幅度降低处理量及处理的复杂程度、有效提高测定的效率和精确度,中心方位角和中心俯仰角的均方根误差仅为角度扩展的十分之一左右,同时避免了多维谱搜索和非线性优化,显著降低了处理的复杂程度。因而,本发明具有方法简单、可靠,可大幅度降低计算处理量及处理的复杂程度,以及测定的效率和精确度高,应用于智能天线系统中、可提高智能天线系统在多径环境中的性能及实际应用的价值等特点。
附图说明
图1是本发明方法的阵列设置示意图;
图2是本发明实施方式仿真运行所得信号一(S1)的中心方位角搜索示轨迹意图;
图3是本发明实施方式仿真运行所得信号一(S1)的中心俯仰角搜索示轨迹意图;
图4是本发明实施方式仿真运行所得信号二(S2)的中心方位角搜索示轨迹意图;
图5是本发明实施方式仿真运行所得信号二(S2)的中心俯仰角搜索示轨迹意图。
具体实施方式
下面结合附图对技术方案的实施作进一步的详细描述。
步骤1.设置立体天线阵列:
设置一个如图1所示的立体天线阵列,它由两个圆阵L1和L2组成,每个圆阵分别设有K=10个阵元,信号波长λ为0.375m,各圆阵半径都为r=λ/(4sin(π/10))=0.3034m;圆阵L1设于x-y平面内,圆阵L2平行于圆阵L1设置;两圆阵间距为d=0.1λ=0.0375m。本实施方式有q=2个窄带非相干分布信号以不同方向入射到此阵列,各阵元上的噪声为加性高斯白噪声,且噪声与信号不相关;
在t时刻,圆阵L1和L2的接收数据为:
其中,为方向矢量。η=2πγ/λ(λ为信号波长,r为园阵半径),γk=2π(k-1)/10(k=1,2,…,10).si(θ,φ,t;μi)为第i个分布信号的角信号密度函数,是对应的角度参数向量,和分别表示中心方位角、方位角扩展、中心俯仰角和俯仰角扩展。
步骤2.建立两个圆阵数据的协方差矩阵:
在工作时,将圆阵L1和L2的接收数据按下式合并为:
一列,并由z(t)的N=400次采样数据{z(1),z(2),…,z(400)}建立如下协方差矩阵
式中,n为采样的序号。
步骤3.确定子矩阵:
对步骤2所得协方差矩阵做特征分解获得特征值和对应的特征向量,并利用最大的6(3q)个特征值所对应的特征向量e1,e2,…,e6作为列建立特征向量矩阵ES={e1,e2,…,e6}。然后再将所得特征向量矩阵(ES)以单个圆阵的阵元数10为准、将其按前、后10行分成两个子矩阵(ES1和ES2),
分别对应于圆阵L1和L2;
步骤4.确定中心俯仰角初值:
将步骤3所得两个子矩阵(ES1和ES2)通过总体最小二乘法建立这两个子矩阵的旋转矩阵(Ψ),具体步骤如下:
(1)建立矩阵ES12=[ES1,ES2].
利用所得旋转矩阵Ψ的特征值ξi(j=1,2,…6)按下式确定各分布信号中心俯仰角的初值
步骤5.建立圆阵L1的协方差矩阵:
由x(t)的N=400次采样数据{x(1),x(2),…,x(400)}确定如下协方差矩阵:
其中,n为采样的序号;
步骤6.确定伪噪声子矩阵:
步骤7.确定二维谱函数:
利用任意圆阵L1的方向矢量、及方向矢量对方位角的偏导数和对俯仰角的偏导数、建立如下广义方向矢量:
步骤8.确定中心方位角:
步骤9.确定中心俯仰角:
将步骤8所得各中心方位角依次代入步骤7所得二维谱函数中,通过如下一维搜索:
即在本实施方式中,信噪比为20dB、采样数为N=400次,对两个分布信号(S1、S2)进行测定。经过测定所得的两个信号的中心方位角和中心俯仰角的值分别为:
为了评估方法的性能,经过200次测定所得的两个信号的中心方位角和中心俯仰角的统计平均值分别为:
相应的中心方位角和中心俯仰角的均方根误差分别为:
其中:两个信号的实际中心方位角和中心俯仰角分别为:
S1:(θ1,φ1)=(20°,65°),
S2:(θ2,φ2)=(30°,75°)。
相应的方位角和俯仰角扩展分别为:
从上可知,按本实施方式,测得的两个信号的中心方位角和中心俯仰角的均方根误差大约为角度扩展的十分之一。而在同样条件下采用传统的点源方法进行测定,测得的均方根误差分别为:
显然,已大于分布信号的角度扩展。
因此,本发明方法相比于点源测定方法大大提高了测定精度。此外,和已有非相干分布源测定方法相比,本发明避免了多维谱搜索和非线性优化,显著降低了处理的复杂程度。
Claims (7)
1.一种非相干分布式信号二维波达角的测定方法,包括;
步骤1.设置立体天线阵列:在三维空间中设立一个由两个均匀圆阵组成的立体天线阵列;每个圆阵的半径及阵元数均相同,其中一个圆阵设于x-y平面内、另一个圆阵则平行于该圆阵设置,以接收各非相干分布信号的数据矢量;
步骤2.建立两个圆阵数据的协方差矩阵:工作时首先将两圆阵中各阵元每次接收的数据分别合并为一列、建立协方差矩阵;
步骤3.确定子矩阵:对步骤2所得协方差矩阵进行特征分解,并利用最大特征值所对应的特征向量建立特征向量矩阵、最大特征值的个数为分布信号数的3倍,然后再将所得特征向量矩阵以阵元数为行数、将其分为前后两个子矩阵;
步骤4.确定中心俯仰角初值:利用总体最小二乘法建立步骤3所得两个子矩阵的旋转矩阵,并利用该旋转矩阵的特征值确定各分布信号中心俯仰角的初值;
步骤5.建立任意一个圆阵的协方差矩阵:根据所选圆阵各次接收的数据建立该圆阵的协方差矩阵;
步骤6.确定伪噪声子矩阵:对步骤5得到的协方差矩阵进行特征分解、得特征值和特征向量,并根据单个圆阵的阵元数与分布信号个数的1.5-2.5倍值得到伪噪声子矩阵的维数,并由该维数个最小特征值所对应的特征向量建立伪噪声子矩阵;
步骤7.确定二维谱函数:利用任意一个圆阵的方向矢量,及该方向矢量对方位角的偏导数和对俯仰角的偏导数、建立广义方向矢量,并利用该广义方向矢量与伪噪声子矩阵之间的正交性确定二维谱函数;
步骤8.确定中心方位角:将步骤4所得各中心俯仰角的初值依次代入步骤7所得二维谱函数中,并用一维搜索方法确定出各初值所对应的中心方位角;
步骤9.确定中心俯仰角:将步骤8所得各中心方位角依次代入步骤7所得二维谱函数中,仍通过一维搜索方法确定出各中心方位角所对应的中心俯仰角。
2.按权利要求1所述非相干分布式信号二维波达角的测定方法,其特征在于所述每个圆阵的阵元数为K=3,4,…20。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201010618502 CN102175989B (zh) | 2010-12-30 | 2010-12-30 | 一种非相干分布式信号二维波达角的测定方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201010618502 CN102175989B (zh) | 2010-12-30 | 2010-12-30 | 一种非相干分布式信号二维波达角的测定方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102175989A CN102175989A (zh) | 2011-09-07 |
CN102175989B true CN102175989B (zh) | 2013-07-17 |
Family
ID=44519192
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN 201010618502 Expired - Fee Related CN102175989B (zh) | 2010-12-30 | 2010-12-30 | 一种非相干分布式信号二维波达角的测定方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102175989B (zh) |
Families Citing this family (13)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103152093B (zh) * | 2011-12-06 | 2017-10-10 | 中兴通讯股份有限公司 | 二维到达角的获取方法和装置 |
CN102608565B (zh) * | 2012-03-23 | 2014-01-08 | 哈尔滨工程大学 | 一种基于均匀圆阵列的波达方向估计方法 |
CN104833947B (zh) * | 2015-04-03 | 2017-10-31 | 西北大学 | 一种任意阵接收对称虚拟变换2d‑doa分离算法 |
CN105510872B (zh) * | 2016-01-22 | 2017-11-17 | 电子科技大学 | 一种适用于大规模mimo系统的二维波达角的测定方法 |
CN108572347A (zh) * | 2017-03-09 | 2018-09-25 | 上海交通大学 | 基于通信信号信道状态响应信息的面阵二维测角方法及系统 |
CN106802402B (zh) * | 2017-03-09 | 2019-04-19 | 西安电子科技大学 | 基于双层平行圆阵列天线的doa估计方法 |
CN107340512B (zh) * | 2017-06-29 | 2020-02-18 | 电子科技大学 | 一种基于子阵划分的近远场混合源被动定位方法 |
CN107271959B (zh) * | 2017-07-05 | 2021-01-26 | 江苏先安科技有限公司 | 基于多天线正交频分复用调制信号的定位方法 |
CN107907853A (zh) * | 2017-11-07 | 2018-04-13 | 中国人民解放军信息工程大学 | 一种基于均匀圆阵差分相位的单分布源doa估计方法 |
CN109188346B (zh) * | 2018-08-31 | 2023-03-10 | 西安电子科技大学 | 大规模均匀圆柱阵列单快拍doa估计方法 |
CN111722177B (zh) * | 2019-03-22 | 2023-06-16 | 成都信息工程大学 | 确定辐射源定向误差的方法 |
CN112051540B (zh) * | 2020-09-11 | 2021-06-01 | 成都大学 | 一种快速高精度测向方法 |
CN112968303B (zh) * | 2021-04-01 | 2022-04-01 | 珠海极海半导体有限公司 | 阵列天线、定位方法、定位系统、ble定位装置及ble设备 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR101351123B1 (ko) * | 2007-03-23 | 2014-01-14 | 삼성전자주식회사 | 가간섭성 소스의 도래방향각 추정을 위한 방법 및 장치 |
CN101908918B (zh) * | 2010-07-26 | 2015-04-01 | 重庆大学 | 无线通信接收机中的一种波束合成方法 |
-
2010
- 2010-12-30 CN CN 201010618502 patent/CN102175989B/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN102175989A (zh) | 2011-09-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102175989B (zh) | 一种非相干分布式信号二维波达角的测定方法 | |
CN106054123B (zh) | 一种稀疏l阵及其二维doa估计方法 | |
CN102169170B (zh) | 一种相干分布式信号二维波达角的测定方法 | |
CN104935366B (zh) | 毫米波通信中的波束搜索方法 | |
CN106019213B (zh) | 一种部分稀疏l阵及其二维doa估计方法 | |
CN106526530B (zh) | 基于传播算子的2-l型阵列二维doa估计算法 | |
CN106972882B (zh) | 基于虚拟域空间功率谱估计的互质阵列自适应波束成形方法 | |
CN104023395B (zh) | 大规模mimo中基于波束空间转换的散射信源定位方法 | |
CN106483493A (zh) | 一种稀疏双平行线阵及二维波达方向估计方法 | |
CN106526531A (zh) | 基于三维天线阵列的改进传播算子二维doa估计算法 | |
US20160156399A1 (en) | Method and device for acquiring channel information | |
CN106291453B (zh) | 一种基于波束空间变换的二维波达角测定方法 | |
CN112929962B (zh) | 定位方法、装置、计算机设备和存储介质 | |
CN107104720A (zh) | 基于协方差矩阵虚拟域离散化重建的互质阵列自适应波束成形方法 | |
CN105933254B (zh) | 毫米波多小区多天线系统中基于波束空间的信道估计方法 | |
CN109782218A (zh) | 一种基于双平行天线阵的相干分布式非圆信号doa估计方法 | |
CN107290732A (zh) | 一种量子大爆炸的单基地mimo雷达测向方法 | |
CN112255629A (zh) | 基于联合uca阵列的序贯esprit二维不相干分布源参数估计方法 | |
CN104506256B (zh) | 一种mimo多天线系统性能评估方法及多天线系统 | |
CN109411890B (zh) | 一种基于多天线的单站高精度三维定位方法 | |
CN108572347A (zh) | 基于通信信号信道状态响应信息的面阵二维测角方法及系统 | |
CN105510872B (zh) | 一种适用于大规模mimo系统的二维波达角的测定方法 | |
Wang et al. | DOA estimation of smart antenna signal based on MUSIC algorithm | |
CN111614386B (zh) | 一种快速搜索方法 | |
CN101686074B (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 | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20130717 Termination date: 20151230 |
|
EXPY | Termination of patent right or utility model |