CN112533284B - 一种基于到达角的近远场统一定位方法 - Google Patents

一种基于到达角的近远场统一定位方法 Download PDF

Info

Publication number
CN112533284B
CN112533284B CN202011020829.XA CN202011020829A CN112533284B CN 112533284 B CN112533284 B CN 112533284B CN 202011020829 A CN202011020829 A CN 202011020829A CN 112533284 B CN112533284 B CN 112533284B
Authority
CN
China
Prior art keywords
target source
column
representing
coordinate system
row
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
CN202011020829.XA
Other languages
English (en)
Other versions
CN112533284A (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.)
Ningbo University
Original Assignee
Ningbo 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 Ningbo University filed Critical Ningbo University
Publication of CN112533284A publication Critical patent/CN112533284A/zh
Application granted granted Critical
Publication of CN112533284B publication Critical patent/CN112533284B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • HELECTRICITY
    • H04ELECTRIC COMMUNICATION TECHNIQUE
    • H04WWIRELESS COMMUNICATION NETWORKS
    • H04W64/00Locating users or terminals or network equipment for network management purposes, e.g. mobility management
    • HELECTRICITY
    • H04ELECTRIC COMMUNICATION TECHNIQUE
    • H04WWIRELESS COMMUNICATION NETWORKS
    • H04W84/00Network topologies
    • H04W84/18Self-organising networks, e.g. ad-hoc networks or sensor networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Signal Processing (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种基于到达角的近远场统一定位方法,其通过建立修正极坐标系来重新表示目标源的位置,建立近远场统一的定位模型;通过泰勒展开将基于到达角的近场非线性定位3D模型转化为线性模型,利用凸松弛技术建立混合半正定/二阶锥规划问题,求解得到近场时目标源的坐标位置估计及远场时目标源的方位角估计和俯仰角估计;通过求得基于到达角的近场非线性定位3D模型的偏差均值的理论值,获得近场时目标源的坐标位置的近似无偏估计值;优点是不需要知道目标源位于近场还是远场的先验知识,能够得到近似无偏估计的目标源位置,且计算量较小。

Description

一种基于到达角的近远场统一定位方法
技术领域
本发明涉及一种近远场统一定位方法,尤其是涉及一种无线传感器网络中基于到达角(Angle-of-Arrive,AOA)的近远场统一定位方法。
背景技术
目标定位是无线传感器网络中一项十分重要的技术。待定位的目标被称为源。在传统的目标定位中,近场定位和远场定位被视为两种不同的定位场景。在近场定位中,目标源距离传感器阵列较近,传感器到目标源的信号传播路径交叉得到目标源位置,因此近场定位通常被视为点定位。在远场定位中,目标源距离传感器阵列较远,传感器到目标源的信号传播路径近似于平行线,无法交叉获得目标源位置,只能估计信号到达方向。近场定位方法往往无法定位远场目标,会产生阈值效应;远场定位方法在目标源不够远时,会出现较大的偏差。目前,业内人员提出了基于到达时间差的近远场统一定位模型(参见G.Wang andK.C.Ho,"Convex Relaxation Methods for Unified Near-Field and Far-Field TDOA-Based Localization"(基于TDOA的统一近远场定位的凸松弛方法).IEEE Transactionson Wireless Communications[J].,vol.18,no.4,pp.2346-2360),但并未研究过基于到达角(Angle-of-Arrive,AOA)的近远场统一定位模型。
在实际应用中,基于到达角的定位由于具有较高的定位精度和不需要严格的时间同步,因而被广泛地应用。但是,基于到达角的定位的测量模型是高度非线性的。对非线性问题求解可以利用尝试所有可能的目标源位置以实现最优估计的直接网格搜索方法,然而这样虽然可以获得最优值,但是计算量庞大。非线性问题的另一种求解方法是通过迭代算法,迭代算法需要一个可靠的初始值开启迭代,初始值的选取决定了整个迭代算法性能的好坏,而合适的初始值的选取是比较困难的。
因此,有必要研究一种定位精度高、计算量小、不需要选取合适的迭代初始值的基于到达角的近远场统一定位方法。
发明内容
本发明所要解决的技术问题是提供一种基于到达角的近远场统一定位方法,其定位精度高,计算量小,无需选取迭代初始值。
本发明解决上述技术问题所采用的技术方案为:一种基于到达角的近远场统一定位方法,其特征在于包括以下步骤:
步骤一:在无线传感器网络中建立一个笛卡尔坐标系作为参考坐标系;设定无线传感器网络中存在一个用于发射测量信号的目标源和N个用于接收测量信号的接收传感器;将N个接收传感器在参考坐标系中的坐标位置的真实值对应记为s1,s2,…,si,...,sN,将目标源在参考坐标系中的坐标位置的真实值记为uo;其中,N>1,s1表示第1个接收传感器在参考坐标系中的坐标位置的真实值,s2表示第2个接收传感器在参考坐标系中的坐标位置的真实值,si表示第i个接收传感器在参考坐标系中的坐标位置的真实值,sN表示第N个接收传感器在参考坐标系中的坐标位置的真实值,si=(xi,yi,zi),xi,yi,zi对应表示si在x轴上的分量、si在y轴上的分量、si在z轴上的分量,1≤i≤N,uo=(xo,yo,zo),xo,yo,zo对应表示uo在x轴上的分量、uo在y轴上的分量、uo在z轴上的分量;
步骤二:将基于到达角的近场非线性定位3D模型描述为:
Figure BDA0002700554240000021
并将第i个接收传感器到目标源的带噪声的方位角测量值记为θi,将第i个接收传感器到目标源的带噪声的俯仰角测量值记为φi
Figure BDA0002700554240000022
然后将
Figure BDA0002700554240000023
代入基于到达角的近场非线性定位3D模型中,得到
Figure BDA0002700554240000031
再分别对
Figure BDA0002700554240000032
中的两个三角函数式用泰勒展开式进行展开,并忽略泰勒展开式中的二阶噪声项使基于到达角的近场非线性定位3D模型转化为线性表达式,描述为:
Figure BDA0002700554240000033
其中,
Figure BDA0002700554240000034
表示第i个接收传感器到目标源的方位角的真实值,
Figure BDA0002700554240000035
表示表示第i个接收传感器到目标源的俯仰角的真实值,
Figure BDA0002700554240000036
表示第i个接收传感器到目标源的方位角的噪声,
Figure BDA0002700554240000037
表示第i个接收传感器到目标源的俯仰角的噪声,
Figure BDA0002700554240000038
Figure BDA0002700554240000039
均服从均值为零和方差为σ2的高斯分布,σ2即为噪声功率;
步骤三:建立一个修正极坐标系,其包括方位角分量、俯仰角分量、长度倒数分量;然后获取目标源在修正极坐标系中的坐标位置的真实值,记为
Figure BDA00027005542400000310
Figure BDA00027005542400000311
其中,θo
Figure BDA00027005542400000312
的方位角分量,θo表示在参考坐标系中目标源相对于原点的方位角的真实值,φo
Figure BDA00027005542400000313
的俯仰角分量,φo表示在参考坐标系中目标源相对于原点的俯仰角的真实值,go
Figure BDA00027005542400000314
的长度倒数分量,go表示在参考坐标系中目标源相对于原点的距离的倒数的真实值,
Figure BDA0002700554240000041
符号“|| ||”为求欧几里德范数符号;
步骤四:将
Figure BDA0002700554240000042
代入线性表达式中,得到线性模型,描述为:
Figure BDA0002700554240000043
步骤五:将线性模型转化为一个约束问题,描述为:
Figure BDA0002700554240000044
其中,min表示最小化,s.t.表示“受约束于……”,v为引入的中间向量,v=[g,ρT,t1,...,ti,...,tN]T,g表示在参考坐标系中目标源相对于原点的距离的倒数的估计值,ρ为引入的中间向量,ρ=[cosθcosφ,sinθcosφ,sinφ],θ表示在参考坐标系中目标源相对于原点的方位角的估计值,φ表示在参考坐标系中目标源相对于原点的俯仰角的估计值,t1,...,ti,...,tN均为引入的中间变量,t1,...,ti,...,tN中的t1和tN均通过ti=||ρ(1:2)-si(1:2)g||计算得到,
Figure BDA0002700554240000045
x1、y1、z1对应表示s1在x轴上的分量、s1在y轴上的分量、s1在z轴上的分量,xN、yN、zN对应表示sN在x轴上的分量、sN在y轴上的分量、sN在z轴上的分量,θ1表示第1个接收传感器到目标源的带噪声的方位角测量值,θN表示第N个接收传感器到目标源的带噪声的方位角测量值,φ1表示第1个接收传感器到目标源的带噪声的俯仰角测量值,φN表示第N个接收传感器到目标源的带噪声的俯仰角测量值,W表示加权矩阵,W=BTQ-1B,B和Q均为引入的中间矩阵,B=diag(t1,...,ti,...,tN),diag(t1,...,ti,...,tN)表示以t1,...,ti,...,tN为对角元素的矩阵,
Figure BDA0002700554240000051
Qθ表示nθ的协方差矩阵,Qφ表示nφ的协方差矩阵,0N×N表示维数为N×N且元素全为0的矩阵,nθ表示由
Figure BDA0002700554240000052
构成的向量,nφ表示由
Figure BDA0002700554240000053
构成的向量,ρ(1:2)表示由ρ中的第1个元素和第2个元素构成的新向量,si(1:2)表示由si中的第1个元素和第2个元素构成的新向量,si(1:2)=[xi,yi],上标“T”表示向量或矩阵的转置;
步骤六:引入一个中间矩阵V,令V=vvT,将约束问题等价转化为一个优化问题,描述为:
Figure BDA0002700554240000054
其中,Tr()表示求矩阵的迹,F、D0、Di均为引入的中间矩阵,F=ATQ-1A,
Figure BDA0002700554240000055
C0为引入的中间矩阵,C0=[03×1,I3,03×N],03×1表示维数为3×1且元素全为0的列向量,I3表示维数为3×3的单位矩阵,03×N表示维数为3×N且元素全为0的矩阵,
Figure BDA0002700554240000056
Ci为引入的中间矩阵,Ci=[-si(1:2),I2,02×N],I2表示维数为2×2的单位矩阵,02×N表示维数为2×N且元素全为0的矩阵,V(4+i,4+i)表示V中第4+i行第4+i列的元素,V≥0表示V是半正定的,rank()表示求矩阵的秩;
步骤七:将优化问题松弛为凸的混合半正定/二阶锥规划问题,描述为:
Figure BDA0002700554240000061
其中,V(2:3,1)表示V的第1列中从第2行到第3行所有元素构成的新向量,V(1,1)表示V中第1行第1列的元素,V(4+i,1)表示V中第4+i行第1列的元素,V(2:3,j+4)表示V的第j+4列中从第2行到第3行所有元素构成的新向量,V(1,j+4)表示V中第1行第j+4列的元素,V(4+i,j+4)表示V中第4+i行第j+4列的元素;
步骤八:采用内点法软件求解混合半正定/二阶锥规划问题,求解得到V的最优解,记为
Figure BDA0002700554240000062
然后对
Figure BDA0002700554240000063
做特征值分解,即
Figure BDA0002700554240000064
再根据
Figure BDA0002700554240000065
得到g、θ、φ各自的最优解,对应记为
Figure BDA0002700554240000066
Figure BDA0002700554240000067
其中,
Figure BDA0002700554240000068
表示v的最优解,
Figure BDA0002700554240000069
对应表示
Figure BDA00027005542400000610
中的第1个元素、第2个元素、第3个元素、第4个元素;
步骤九:计算在近场情况下目标源在参考坐标系中的坐标位置的估计,记为
Figure BDA00027005542400000611
Figure BDA00027005542400000612
其中,
Figure BDA00027005542400000613
在远场情况下,
Figure BDA00027005542400000614
为在参考坐标系中目标源相对于原点的方位角的估计,
Figure BDA00027005542400000615
为在参考坐标系中目标源相对于原点的俯仰角的估计,
Figure BDA00027005542400000616
没有意义;
步骤十:在步骤二中的泰勒展开式中保留其二阶噪声项,然后采用拉格朗日乘子法对泰勒展开式进行求解,求得基于到达角的近场非线性定位3D模型的偏差均值的理论值,记为
Figure BDA00027005542400000617
Figure BDA00027005542400000618
其中,
Figure BDA0002700554240000071
表示基于到达角的近场非线性定位3D模型的偏差,
Figure BDA0002700554240000072
Go为引入的中间矩阵,Go=(Uo)TWUo,Uo为U不带误差的矩阵,Uo=U-ΔU,
Figure BDA0002700554240000073
Figure BDA0002700554240000074
表示目标源在修正极坐标系中的坐标位置的变量,ΔU为误差矩阵,ΔU的第1列元素ΔU(:,1)为
Figure BDA0002700554240000075
ΔU的第2列元素ΔU(:,2)为
Figure BDA0002700554240000076
ΔU的第3列元素ΔU(:,3)为
Figure BDA0002700554240000077
θ1 o表示第1个接收传感器到目标源的方位角的真实值,φ1 o表示第1个接收传感器到目标源的俯仰角的真实值,p1、pi
Figure BDA0002700554240000078
f1、fi
Figure BDA0002700554240000079
均为引入的中间变量,pi=xigo sinθo cosφo-yigo cosθo cosφo,p1根据pi的计算公式计算得到,
Figure BDA0002700554240000081
Figure BDA0002700554240000082
根据
Figure BDA0002700554240000083
的计算公式计算得到,
Figure BDA0002700554240000084
f1根据fi的计算公式计算得到,
Figure BDA0002700554240000085
ρo(1:2)表示由ρ的真实值ρo中的第1个元素和第2个元素构成的新向量,
Figure BDA0002700554240000086
根据
Figure BDA0002700554240000087
的计算公式计算得到,E[n⊙n]=[01×N,Q],01×N表示维数为1×N且元素全为0的向量,符号“⊙”表示两个向量对应元素相乘,n=[nθ,nφ],L为引入的中间矩阵,
Figure BDA0002700554240000088
Figure BDA0002700554240000089
表示以
Figure BDA00027005542400000810
为对角元素的矩阵,
Figure BDA00027005542400000811
Figure BDA00027005542400000812
表示
Figure BDA00027005542400000813
的均值的理论值,
Figure BDA00027005542400000814
Figure BDA00027005542400000815
表示
Figure BDA00027005542400000816
的均值的理论值,
Figure BDA00027005542400000817
Figure BDA00027005542400000818
表示
Figure BDA00027005542400000819
的第d列,
Figure BDA00027005542400000820
表示
Figure BDA00027005542400000821
的第d+N列,qd表示Q的第d列,qd+N表示Q的第d+N列,
Figure BDA00027005542400000822
表示
Figure BDA00027005542400000823
的第k列,ΔU(d,1)表示ΔU中第d行第1列的元素,ΔU(d,4)表示ΔU中第d行第4列的元素,ΔU(d,1+N)表示ΔU中第d行第1+N列的元素,ΔU(d,4+N)表示ΔU中第d行第4+N列的元素,ΔU(1,1:N)表示由ΔU的第1行中第1列到第N列的所有元素构成的新向量,ΔU(N,1:N)表示由ΔU的第N行中第1列到第N列的所有元素构成的新向量,ΔU(1,1+N:2N)表示由ΔU的第1行中第1+N列到第2N列的所有元素构成的新向量,ΔU(N,1+N:2N)表示由ΔU的第N行中第1+N列到第2N列的所有元素构成的新向量,Q(1,k)表示Q中第1行第k列的元素,Q(N,k)表示Q中第N行第k列的元素,Q(N+1,k)表示Q中第N+1行第k列的元素,Q(2N,k)表示Q中第2N行第k列的元素;
步骤十一:将在近场情况下目标源在参考坐标系中的坐标位置的估计减去基于到达角的近场非线性定位3D模型的偏差均值的理论值,得到目标源的坐标位置的近似无偏估计值,记为
Figure BDA0002700554240000091
Figure BDA0002700554240000092
与现有技术相比,本发明的优点在于:
本发明方法通过泰勒展开处理基于到达角的近场非线性定位3D模型,转化为线性表达式,最终构建凸的混合半正定/二阶锥规划问题来求解目标源的坐标位置的估计值,由于求解凸的混合半正定/二阶锥规划问题可以采用内点法,因此无需选取可靠的初始值;本发明方法相较于传统的全局搜索方法降低了计算量;本发明方法能够在不需要目标源位于近场还是远场的先验知识,就能实现近远场统一定位,而且还提供了一个近似无偏的估计值,定位精度高。
附图说明
图1为本发明方法的总体实现框图;
图2a为本发明方法(SDP-MPR-BR)与现有的AOA定位方法(BRPLE)在噪声功率σ2=0.001的情况下目标源位置的角度对数均方误差(MSE)随在参考坐标系下目标源到原点的距离变化的示意图;
图2b为本发明方法(SDP-MPR-BR)与现有的AOA定位方法(BRPLE)在噪声功率σ2=0.001的情况下目标源位置的反距离对数均方误差(MSE)随在参考坐标系下目标源到原点的距离变化的示意图;
图3a为本发明方法(SDP-MPR-BR)与现有的AOA定位方法(BRPLE)在噪声功率σ2=0.001的情况下目标源位置的角度对数偏差(Bias)随在参考坐标系下目标源到原点的距离变化的示意图;
图3b为本发明方法(SDP-MPR-BR)与现有的AOA定位方法(BRPLE)在噪声功率σ2=0.001的情况下目标源位置的反距离对数偏差(Bias)随在参考坐标系下目标源到原点的距离变化的示意图。
具体实施方式
以下结合附图实施例对本发明作进一步详细描述。
本发明提出的一种基于到达角的近远场统一定位方法,其总体实现框图如图1所示,其包括以下步骤:
步骤一:在无线传感器网络中建立一个笛卡尔坐标系作为参考坐标系;设定无线传感器网络中存在一个用于发射测量信号的目标源和N个用于接收测量信号的接收传感器;将N个接收传感器在参考坐标系中的坐标位置的真实值对应记为s1,s2,…,si,…,sN,将目标源在参考坐标系中的坐标位置的真实值记为uo;其中,N>1,在实验中取N=6,s1表示第1个接收传感器在参考坐标系中的坐标位置的真实值,s2表示第2个接收传感器在参考坐标系中的坐标位置的真实值,si表示第i个接收传感器在参考坐标系中的坐标位置的真实值,sN表示第N个接收传感器在参考坐标系中的坐标位置的真实值,si=(xi,yi,zi),xi,yi,zi对应表示si在x轴上的分量、si在y轴上的分量、si在z轴上的分量,1≤i≤N,uo=(xo,yo,zo),xo,yo,zo对应表示uo在x轴上的分量、uo在y轴上的分量、uo在z轴上的分量。
步骤二:将已知的基于到达角的近场非线性定位3D模型描述为:
Figure BDA0002700554240000111
并将第i个接收传感器到目标源的带噪声的方位角测量值记为θi,将第i个接收传感器到目标源的带噪声的俯仰角测量值记为φi
Figure BDA0002700554240000112
然后将
Figure BDA0002700554240000113
代入基于到达角的近场非线性定位3D模型中,得到
Figure BDA0002700554240000114
再分别对
Figure BDA0002700554240000115
中的两个三角函数式用泰勒展开式进行展开,并忽略泰勒展开式中的二阶噪声项使基于到达角的近场非线性定位3D模型转化为线性表达式,描述为:
Figure BDA0002700554240000116
其中
Figure BDA0002700554240000117
表示第i个接收传感器到目标源的方位角的真实值,
Figure BDA0002700554240000118
表示表示第i个接收传感器到目标源的俯仰角的真实值,
Figure BDA0002700554240000119
表示第i个接收传感器到目标源的方位角的噪声,
Figure BDA00027005542400001110
表示第i个接收传感器到目标源的俯仰角的噪声,
Figure BDA00027005542400001111
Figure BDA00027005542400001112
均服从均值为零和方差为σ2的高斯分布,σ2根据环境自行设定,σ2即为噪声功率。
步骤三:建立一个修正极坐标系,其包括方位角分量、俯仰角分量、长度倒数分量;然后获取目标源在修正极坐标系中的坐标位置的真实值,记为
Figure BDA0002700554240000121
Figure BDA0002700554240000122
其中,θo
Figure BDA0002700554240000123
的方位角分量,θo表示在参考坐标系中目标源相对于原点的方位角的真实值,φo
Figure BDA0002700554240000124
的俯仰角分量,φo表示在参考坐标系中目标源相对于原点的俯仰角的真实值,go
Figure BDA0002700554240000125
的长度倒数分量,go表示在参考坐标系中目标源相对于原点的距离的倒数的真实值,
Figure BDA0002700554240000126
符号“|| |”为求欧几里德范数符号。
步骤四:将
Figure BDA0002700554240000127
代入线性表达式中,得到线性模型,描述为:
Figure BDA0002700554240000128
步骤五:将线性模型转化为一个约束问题,描述为:
Figure BDA0002700554240000129
其中,min表示最小化,s.t.表示“受约束于……”,v为引入的中间向量,v=[g,ρT,t1,…,ti,…,tN]T,g表示在参考坐标系中目标源相对于原点的距离的倒数的估计值,ρ为引入的中间向量,ρ=[cosθcosφ,sinθcosφ,sinφ],θ表示在参考坐标系中目标源相对于原点的方位角的估计值,φ表示在参考坐标系中目标源相对于原点的俯仰角的估计值,t1,…,ti,…,tN均为引入的中间变量,t1,…,ti,…,tN中的t1和tN均通过ti=||ρ(1:2)-si(1:2)g||计算得到,
Figure BDA0002700554240000131
x1、y1、z1对应表示s1在x轴上的分量、s1在y轴上的分量、s1在z轴上的分量,xN、yN、zN对应表示sN在x轴上的分量、sN在y轴上的分量、sN在z轴上的分量,θ1表示第1个接收传感器到目标源的带噪声的方位角测量值,θN表示第N个接收传感器到目标源的带噪声的方位角测量值,φ1表示第1个接收传感器到目标源的带噪声的俯仰角测量值,φN表示第N个接收传感器到目标源的带噪声的俯仰角测量值,W表示加权矩阵,W=BTQ-1B,B和Q均为引入的中间矩阵,B=diag(t1,...,ti,...,tN),diag(t1,...,ti,...,tN)表示以t1,...,ti,...,tN为对角元素的矩阵,
Figure BDA0002700554240000132
Qθ表示nθ的协方差矩阵,Qφ表示nφ的协方差矩阵,0N×N表示维数为N×N且元素全为0的矩阵,nθ表示由
Figure BDA0002700554240000133
构成的向量,nφ表示由
Figure BDA0002700554240000134
构成的向量,ρ(1:2)表示由ρ中的第1个元素和第2个元素构成的新向量,si(1:2)表示由si中的第1个元素和第2个元素构成的新向量,si(1:2)=[xi,yi],上标“T”表示向量或矩阵的转置。
步骤六:引入一个中间矩阵V,令V=vvT,将约束问题等价转化为一个优化问题,描述为:
Figure BDA0002700554240000135
其中,Tr()表示求矩阵的迹,F、D0、Di均为引入的中间矩阵,F=ATQ-1A,
Figure BDA0002700554240000136
C0为引入的中间矩阵,C0=[03×1,I3,03×N],03×1表示维数为3×1且元素全为0的列向量,I3表示维数为3×3的单位矩阵,03×N表示维数为3×N且元素全为0的矩阵,
Figure BDA0002700554240000141
Ci为引入的中间矩阵,Ci=[-si(1:2),I2,02×N],I2表示维数为2×2的单位矩阵,02×N表示维数为2×N且元素全为0的矩阵,V(4+i,4+i)表示V中第4+i行第4+i列的元素,V≥0表示V是半正定的,rank()表示求矩阵的秩。
步骤七:将优化问题松弛为凸的混合半正定/二阶锥规划问题,描述为:
Figure BDA0002700554240000142
其中,V(2:3,1)表示V的第1列中从第2行到第3行所有元素构成的新向量,V(1,1)表示V中第1行第1列的元素,V(4+i,1)表示V中第4+i行第1列的元素,V(2:3,j+4)表示V的第j+4列中从第2行到第3行所有元素构成的新向量,V(1,j+4)表示V中第1行第j+4列的元素,V(4+i,j+4)表示V中第4+i行第j+4列的元素。
步骤八:采用常见的内点法软件求解混合半正定/二阶锥规划问题,求解得到V的最优解,记为
Figure BDA0002700554240000143
然后对
Figure BDA0002700554240000144
做特征值分解,即
Figure BDA0002700554240000145
再根据
Figure BDA0002700554240000146
得到g、θ、φ各自的最优解,对应记为
Figure BDA0002700554240000147
Figure BDA0002700554240000148
其中,
Figure BDA0002700554240000149
表示v的最优解,
Figure BDA00027005542400001410
对应表示
Figure BDA00027005542400001411
中的第1个元素、第2个元素、第3个元素、第4个元素。
步骤九:计算在近场情况下目标源在参考坐标系中的坐标位置的估计,记为
Figure BDA00027005542400001412
Figure BDA00027005542400001413
其中,
Figure BDA00027005542400001414
在远场情况下,
Figure BDA00027005542400001415
为在参考坐标系中目标源相对于原点的方位角的估计,
Figure BDA0002700554240000151
为在参考坐标系中目标源相对于原点的俯仰角的估计,
Figure BDA0002700554240000152
没有意义。
步骤十:在步骤二中的泰勒展开式中保留其二阶噪声项,然后采用拉格朗日乘子法对泰勒展开式进行求解,求得基于到达角的近场非线性定位3D模型的偏差均值的理论值,记为
Figure BDA0002700554240000153
Figure BDA0002700554240000154
其中,
Figure BDA0002700554240000155
表示基于到达角的近场非线性定位3D模型的偏差,
Figure BDA0002700554240000156
Go为引入的中间矩阵,Go=(Uo)TWUo,Uo为U不带误差的矩阵,Uo=U-ΔU,
Figure BDA0002700554240000157
Figure BDA00027005542400001510
表示目标源在修正极坐标系中的坐标位置的变量,ΔU为误差矩阵,ΔU的第1列元素ΔU(:,1)为
Figure BDA0002700554240000158
ΔU的第2列元素ΔU(:,2)为
Figure BDA0002700554240000159
ΔU的第3列元素ΔU(:,3)为
Figure BDA0002700554240000161
θ1 o表示第1个接收传感器到目标源的方位角的真实值,φ1 o表示第1个接收传感器到目标源的俯仰角的真实值,p1、pi
Figure BDA0002700554240000162
f1、fi
Figure BDA0002700554240000163
均为引入的中间变量,pi=xigo sinθo cosφo-yigo cosθo cosφo,p1根据pi的计算公式计算得到,
Figure BDA0002700554240000164
Figure BDA0002700554240000165
根据
Figure BDA0002700554240000166
的计算公式计算得到,
Figure BDA0002700554240000167
f1根据fi的计算公式计算得到,
Figure BDA0002700554240000168
ρo(1:2)表示由ρ的真实值ρo中的第1个元素和第2个元素构成的新向量,
Figure BDA0002700554240000169
根据
Figure BDA00027005542400001610
的计算公式计算得到,E[n⊙n]=[01×N,Q],01×N表示维数为1×N且元素全为0的向量,符号“⊙”表示两个向量对应元素相乘,n=[nθ,nφ],L为引入的中间矩阵,
Figure BDA00027005542400001611
Figure BDA00027005542400001612
表示以
Figure BDA00027005542400001613
为对角元素的矩阵,
Figure BDA00027005542400001614
Figure BDA00027005542400001615
表示
Figure BDA00027005542400001616
的均值的理论值,
Figure BDA00027005542400001617
Figure BDA00027005542400001618
表示
Figure BDA0002700554240000171
的均值的理论值,
Figure BDA0002700554240000172
Figure BDA0002700554240000173
表示
Figure BDA0002700554240000174
的第d列,
Figure BDA0002700554240000175
表示
Figure BDA0002700554240000176
的第d+N列,qd表示Q的第d列,qd+N表示Q的第d+N列,
Figure BDA0002700554240000177
表示
Figure BDA0002700554240000178
的第k列,ΔU(d,1)表示ΔU中第d行第1列的元素,ΔU(d,4)表示ΔU中第d行第4列的元素,ΔU(d,1+N)表示ΔU中第d行第1+N列的元素,ΔU(d,4+N)表示ΔU中第d行第4+N列的元素,ΔU(1,1:N)表示由ΔU的第1行中第1列到第N列的所有元素构成的新向量,ΔU(N,1:N)表示由ΔU的第N行中第1列到第N列的所有元素构成的新向量,ΔU(1,1+N:2N)表示由ΔU的第1行中第1+N列到第2N列的所有元素构成的新向量,ΔU(N,1+N:2N)表示由ΔU的第N行中第1+N列到第2N列的所有元素构成的新向量,Q(1,k)表示Q中第1行第k列的元素,Q(N,k)表示Q中第N行第k列的元素,Q(N+1,k)表示Q中第N+1行第k列的元素,Q(2N,k)表示Q中第2N行第k列的元素。
步骤十一:将在近场情况下目标源在参考坐标系中的坐标位置的估计减去基于到达角的近场非线性定位3D模型的偏差均值的理论值,得到目标源的坐标位置的近似无偏估计值,记为
Figure BDA0002700554240000179
Figure BDA00027005542400001710
为验证本发明方法的可行性和有效性,对本发明方法进行仿真试验。
假设无线传感器网络中有6个接收传感器,其在参考坐标系中的坐标位置的真实值如表1所列。在参考坐标系中目标源相对于原点的方位角的真实值设为θo=35.12°,在参考坐标系中目标源相对于原点的俯仰角的真实值设为φo=-62.95°。所有接收传感器到目标源的方位角的噪声的协方差矩阵Qθ=σ2IN,所有接收传感器到目标源的俯仰角的噪声的协方差矩阵Qφ=σ2IN,其中,σ2为噪声功率,IN为维数为N×N的单位矩阵。
表1接收传感器的位置
Figure BDA0002700554240000181
测试本发明方法在固定噪声功率的情况下,定位性能随在参考坐标系下目标源到原点的距离变化的情况。用于对比的方法为Y.Wang and K.C.Ho.An asymptoticallyefficient estimator in closed-form for 3-D AOA localization using a sensornetwork(传感器网络中的3-D AOA定位的闭式渐近无偏估计方法).IEEE Transactions onWireless Communications[J].2015,14(12):6524-6535.中公开的偏差减小的伪线性估计器(Bias Reduction Pseudolinear Estimator,BRPLE)。克拉美罗下界(Cramer-Rao lowerbound,CRLB)为模型理论均方误差下界。
图2a给出了本发明方法(SDP-MPR-BR)与现有的AOA定位方法(BRPLE)在噪声功率σ2=0.001的情况下目标源位置的角度对数均方误差(MSE)随在参考坐标系下目标源到原点的距离变化的示意图;图2b给出了本发明方法(SDP-MPR-BR)与现有的AOA定位方法(BRPLE)在噪声功率σ2=0.001的情况下目标源位置的反距离对数均方误差(MSE)随在参考坐标系下目标源到原点的距离变化的示意图。图3a给出了本发明方法(SDP-MPR-BR)与现有的AOA定位方法(BRPLE)在噪声功率σ2=0.001的情况下目标源位置的角度对数偏差(Bias)随在参考坐标系下目标源到原点的距离变化的示意图;图3b给出了本发明方法(SDP-MPR-BR)与现有的AOA定位方法(BRPLE)在噪声功率σ2=0.001的情况下目标源位置的反距离对数偏差(Bias)随在参考坐标系下目标源到原点的距离变化的示意图。从图2a、图2b、图3a和图3b中可以看出,近远场的分界点在150左右,现有的AOA定位方法(BRPLE)仅在近场时定位性能良好,但无法估计目标源处于远场时的DOA,相比之下,本发明方法在近场和远场时都有良好的定位性能。
从上述仿真结果可以看出,本发明方法具有良好的定位性能,能够很好地满足高定位精度的需求。

Claims (1)

1.一种基于到达角的近远场统一定位方法,其特征在于包括以下步骤:
步骤一:在无线传感器网络中建立一个笛卡尔坐标系作为参考坐标系;设定无线传感器网络中存在一个用于发射测量信号的目标源和N个用于接收测量信号的接收传感器;将N个接收传感器在参考坐标系中的坐标位置的真实值对应记为s1,s2,...,si,...,sN,将目标源在参考坐标系中的坐标位置的真实值记为uo;其中,N>1,s1表示第1个接收传感器在参考坐标系中的坐标位置的真实值,s2表示第2个接收传感器在参考坐标系中的坐标位置的真实值,si表示第i个接收传感器在参考坐标系中的坐标位置的真实值,sN表示第N个接收传感器在参考坐标系中的坐标位置的真实值,si=(xi,yi,zi),xi,yi,zi对应表示si在x轴上的分量、si在y轴上的分量、si在z轴上的分量,1≤i≤N,uo=(xo,yo,zo),xo,yo,zo对应表示uo在x轴上的分量、uo在y轴上的分量、uo在z轴上的分量;
步骤二:将基于到达角的近场非线性定位3D模型描述为:
Figure FDA0002700554230000011
并将第i个接收传感器到目标源的带噪声的方位角测量值记为θi,将第i个接收传感器到目标源的带噪声的俯仰角测量值记为φi
Figure FDA0002700554230000012
然后将
Figure FDA0002700554230000013
代入基于到达角的近场非线性定位3D模型中,得到
Figure FDA0002700554230000014
再分别对
Figure FDA0002700554230000021
中的两个三角函数式用泰勒展开式进行展开,并忽略泰勒展开式中的二阶噪声项使基于到达角的近场非线性定位3D模型转化为线性表达式,描述为:
Figure FDA0002700554230000022
其中,
Figure FDA0002700554230000023
表示第i个接收传感器到目标源的方位角的真实值,
Figure FDA0002700554230000024
表示第i个接收传感器到目标源的俯仰角的真实值,
Figure FDA0002700554230000025
表示第i个接收传感器到目标源的方位角的噪声,
Figure FDA0002700554230000026
表示第i个接收传感器到目标源的俯仰角的噪声,
Figure FDA0002700554230000027
Figure FDA0002700554230000028
均服从均值为零和方差为σ2的高斯分布,σ2即为噪声功率;
步骤三:建立一个修正极坐标系,其包括方位角分量、俯仰角分量、长度倒数分量;然后获取目标源在修正极坐标系中的坐标位置的真实值,记为
Figure FDA0002700554230000029
其中,θo
Figure FDA00027005542300000210
的方位角分量,θo表示在参考坐标系中目标源相对于原点的方位角的真实值,φo
Figure FDA00027005542300000211
的俯仰角分量,φo表示在参考坐标系中目标源相对于原点的俯仰角的真实值,go
Figure FDA00027005542300000212
的长度倒数分量,go表示在参考坐标系中目标源相对于原点的距离的倒数的真实值,
Figure FDA00027005542300000213
符号“|| ||”为求欧几里德范数符号;
步骤四:将
Figure FDA00027005542300000214
代入线性表达式中,得到线性模型,描述为:
Figure FDA0002700554230000031
步骤五:将线性模型转化为一个约束问题,描述为:
Figure FDA0002700554230000032
其中,min表示最小化,s.t.表示“受约束于……”,v为引入的中间向量,v=[g,ρT,t1,...,ti,...,tN]T,g表示在参考坐标系中目标源相对于原点的距离的倒数的估计值,ρ为引入的中间向量,ρ=[cosθcosφ,sinθcosφ,sinφ],θ表示在参考坐标系中目标源相对于原点的方位角的估计值,φ表示在参考坐标系中目标源相对于原点的俯仰角的估计值,t1,...,ti,...,tN均为引入的中间变量,t1,...,ti,...,tN中的t1和tN均通过ti=||ρ(1:2)-si(1:2)g||计算得到,
Figure FDA0002700554230000033
x1、y1、z1对应表示s1在x轴上的分量、s1在y轴上的分量、s1在z轴上的分量,xN、yN、zN对应表示sN在x轴上的分量、sN在y轴上的分量、sN在z轴上的分量,θ1表示第1个接收传感器到目标源的带噪声的方位角测量值,θN表示第N个接收传感器到目标源的带噪声的方位角测量值,φ1表示第1个接收传感器到目标源的带噪声的俯仰角测量值,φN表示第N个接收传感器到目标源的带噪声的俯仰角测量值,W表示加权矩阵,W=BTQ-1B,B和Q均为引入的中间矩阵,B=diag(t1,...,ti,...,tN),diag(t1,...,ti,...,tN) 表示以t1,...,ti,...,tN为对角元素的矩阵,
Figure FDA0002700554230000041
Qθ表示nθ的协方差矩阵,Qφ表示nφ的协方差矩阵,0N×N表示维数为N×N且元素全为0的矩阵,nθ表示由
Figure FDA0002700554230000042
构成的向量,nφ表示由
Figure FDA0002700554230000043
构成的向量,ρ(1:2)表示由ρ中的第1个元素和第2个元素构成的新向量,si(1:2)表示由si中的第1个元素和第2个元素构成的新向量,si(1:2)=[xi,yi],上标“T”表示向量或矩阵的转置;
步骤六:引入一个中间矩阵V,令V=vvT,将约束问题等价转化为一个优化问题,描述为:
Figure FDA0002700554230000044
其中,Tr()表示求矩阵的迹,F、D0、Di均为引入的中间矩阵,F=ATQ-1A,
Figure FDA0002700554230000045
C0为引入的中间矩阵,C0=[03×1,I3,03×N],03×1表示维数为3×1且元素全为0的列向量,I3表示维数为3×3的单位矩阵,03×N表示维数为3×N且元素全为0的矩阵,
Figure FDA0002700554230000046
Ci为引入的中间矩阵,Ci=[-si(1:2),I2,02×N],I2表示维数为2×2的单位矩阵,02×N表示维数为2×N且元素全为0的矩阵,V(4+i,4+i)表示V中第4+i行第4+i列的元素,V≥0表示V是半正定的,rank()表示求矩阵的秩;
步骤七:将优化问题松弛为凸的混合半正定/二阶锥规划问题,描述为:
Figure FDA0002700554230000047
其中,V(2:3,1)表示V的第1列中从第2行到第3行所有元素构成的新向量,V(1,1)表示V中第1行第1列的元素,V(4+i,1)表示V中第4+i行第1列的元素,V(2:3,j+4)表示V的第j+4列中从第2行到第3行所有元素构成的新向量,V(1,j+4)表示V中第1行第j+4列的元素,V(4+i,j+4)表示V中第4+i行第j+4列的元素;
步骤八:采用内点法软件求解混合半正定/二阶锥规划问题,求解得到V的最优解,记为
Figure FDA0002700554230000051
然后对
Figure FDA0002700554230000052
做特征值分解,即
Figure FDA0002700554230000053
再根据
Figure FDA0002700554230000054
得到g、θ、φ各自的最优解,对应记为
Figure FDA0002700554230000055
其中,
Figure FDA0002700554230000056
表示v的最优解,
Figure FDA0002700554230000057
对应表示
Figure FDA0002700554230000058
中的第1个元素、第2个元素、第3个元素、第4个元素;
步骤九:计算在近场情况下目标源在参考坐标系中的坐标位置的估计,记为
Figure FDA0002700554230000059
Figure FDA00027005542300000510
其中,
Figure FDA00027005542300000511
在远场情况下,
Figure FDA00027005542300000512
为在参考坐标系中目标源相对于原点的方位角的估计,
Figure FDA00027005542300000513
为在参考坐标系中目标源相对于原点的俯仰角的估计,
Figure FDA00027005542300000514
没有意义;
步骤十:在步骤二中的泰勒展开式中保留其二阶噪声项,然后采用拉格朗日乘子法对泰勒展开式进行求解,求得基于到达角的近场非线性定位3D模型的偏差均值的理论值,记为
Figure FDA00027005542300000515
Figure FDA00027005542300000516
其中,
Figure FDA00027005542300000517
表示基于到达角的近场非线性定位3D模型的偏差,
Figure FDA00027005542300000518
Go为引入的中间矩阵,Go=(Uo)TWUo,Uo为U不带误差的矩阵,
Figure FDA00027005542300000519
Figure FDA00027005542300000520
表示目标源在修正极坐标系中的坐标位置的变量,ΔU为误差矩阵,ΔU的第1列元素ΔU(:,1)为
Figure FDA0002700554230000061
ΔU的第2列元素ΔU(:,2)为
Figure FDA0002700554230000062
ΔU的第3列元素ΔU(:,3)为
Figure FDA0002700554230000063
θ1 o表示第1个接收传感器到目标源的方位角的真实值,φ1 o表示第1个接收传感器到目标源的俯仰角的真实值,p1、pi
Figure FDA0002700554230000064
f1、fi
Figure FDA0002700554230000065
均为引入的中间变量,pi=xigosinθocosφo-yigocosθocosφo,p1根据pi的计算公式计算得到,
Figure FDA0002700554230000066
Figure FDA0002700554230000067
根据
Figure FDA0002700554230000068
的计算公式计算得到,
Figure FDA0002700554230000069
f1根据fi的计算公式计算得到,
Figure FDA00027005542300000610
ρo(1:2)表示由ρ的真实值ρo中的第1个元素和第2个元素构成的新向量,
Figure FDA0002700554230000071
根据
Figure FDA0002700554230000072
的计算公式计算得到,E[n⊙n]=[01×N,Q],01×N表示维数为1×N且元素全为0的向量,符号“⊙”表示两个向量对应元素相乘,n=[nθ,nφ],L为引入的中间矩阵,
Figure FDA0002700554230000073
Figure FDA0002700554230000074
表示以01×N,
Figure FDA0002700554230000075
为对角元素的矩阵,
Figure FDA0002700554230000076
Figure FDA0002700554230000077
表示
Figure FDA0002700554230000078
的均值的理论值,
Figure FDA0002700554230000079
Figure FDA00027005542300000710
表示
Figure FDA00027005542300000711
的均值的理论值,
Figure FDA00027005542300000712
Figure FDA00027005542300000713
表示
Figure FDA00027005542300000714
的第d列,
Figure FDA00027005542300000715
表示
Figure FDA00027005542300000716
的第d+N列,qd表示Q的第d列,qd+N表示Q的第d+N列,
Figure FDA00027005542300000717
表示
Figure FDA00027005542300000718
的第k列,ΔU(d,1)表示ΔU中第d行第1列的元素,ΔU(d,4)表示ΔU中第d行第4列的元素,ΔU(d,1+N)表示ΔU中第d行第1+N列的元素,ΔU(d,4+N)表示ΔU中第d行第4+N列的元素,ΔU(1,1:N)表示由ΔU的第1行中第1列到第N列的所有元素构成的新向量,ΔU(N,1:N)表示由ΔU的第N行中第1列到第N列的所有元素构成的新向量,ΔU(1,1+N:2N)表示由ΔU的第1行中第1+N列到第2N列的所有元素构成的新向量,ΔU(N,1+N:2N)表示由ΔU的第N行中第1+N列到第2N列的所有元素构成的新向量,Q(1,k)表示Q中第1行第k列的元素,Q(N,k)表示Q中第N行第k列的元素,Q(N+1,k)表示Q中第N+1行第k列的元素,Q(2N,k)表示Q中第2N行第k列的元素;
步骤十一:将在近场情况下目标源在参考坐标系中的坐标位置的估计减去基于到达角的近场非线性定位3D模型的偏差均值的理论值,得到目标源的坐标位置的近似无偏估计值,记为
Figure FDA0002700554230000081
CN202011020829.XA 2020-09-23 2020-09-25 一种基于到达角的近远场统一定位方法 Active CN112533284B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN2020110081833 2020-09-23
CN202011008183 2020-09-23

Publications (2)

Publication Number Publication Date
CN112533284A CN112533284A (zh) 2021-03-19
CN112533284B true CN112533284B (zh) 2022-05-27

Family

ID=74980339

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011020829.XA Active CN112533284B (zh) 2020-09-23 2020-09-25 一种基于到达角的近远场统一定位方法

Country Status (1)

Country Link
CN (1) CN112533284B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113518358B (zh) * 2021-04-25 2022-05-31 长春理工大学 一种无线传感器网络布测位置误差的校正方法和装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109975745A (zh) * 2019-02-28 2019-07-05 宁波大学 一种基于到达时间差的近远场统一定位方法
CN110095753A (zh) * 2019-05-14 2019-08-06 北京邮电大学 一种基于到达角度aoa测距的定位方法及装置
CN110662163A (zh) * 2019-08-23 2020-01-07 宁波大学 基于rss和aoa的三维无线传感网络协作定位方法
CN110658490A (zh) * 2019-08-23 2020-01-07 宁波大学 基于rss和aoa的三维无线传感网络非协作定位方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7460976B2 (en) * 2004-06-09 2008-12-02 The Board Of Trustees Of The Leland Stanford Jr. University Semi-definite programming method for ad hoc network node localization

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109975745A (zh) * 2019-02-28 2019-07-05 宁波大学 一种基于到达时间差的近远场统一定位方法
CN110095753A (zh) * 2019-05-14 2019-08-06 北京邮电大学 一种基于到达角度aoa测距的定位方法及装置
CN110662163A (zh) * 2019-08-23 2020-01-07 宁波大学 基于rss和aoa的三维无线传感网络协作定位方法
CN110658490A (zh) * 2019-08-23 2020-01-07 宁波大学 基于rss和aoa的三维无线传感网络非协作定位方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
一种基于TDOA/AOA的混合三维定位算法;杨浩等;《南京邮电大学学报(自然科学版)》;20121215(第06期);全文 *
一种基于接收信号强度和到达角的三维组合定位方法;张瑞玲;《舰船电子对抗》;20190425(第02期);全文 *
基于约束总体最小二乘的泰勒级数定位算法;陆剑锋等;《计算机应用与软件》;20191212(第12期);全文 *

Also Published As

Publication number Publication date
CN112533284A (zh) 2021-03-19

Similar Documents

Publication Publication Date Title
CN103969622B (zh) 一种基于多运动接收站的时差定位方法
CN109143152B (zh) 基于张量建模的极化阵列波达方向和极化参数估计方法
CN108668358B (zh) 一种应用于无线传感网络的基于到达时间的协作定位方法
CN109917333B (zh) 融合aoa观测量与tdoa观测量的无源定位方法
CN106405533B (zh) 基于约束加权最小二乘的雷达目标联合同步与定位方法
CN110673089B (zh) 未知视距和非视距分布情况下基于到达时间的定位方法
CN104020440B (zh) 基于l型干涉式线性阵列的二维波达角估计方法
CN113835063B (zh) 一种无人机阵列幅相误差与信号doa联合估计方法
CN104793177B (zh) 基于最小二乘法的麦克风阵列测向方法
CN107703482A (zh) 一种闭式解与迭代算法相结合的aoa定位方法
CN112533284B (zh) 一种基于到达角的近远场统一定位方法
CN111830465B (zh) 二维牛顿正交匹配追踪压缩波束形成方法
CN111273302B (zh) 一种浅海匀速运动目标初始状态估计方法
CN103323810A (zh) 一种l阵方位角和俯仰角配对的信号处理方法
CN113030847A (zh) 一种用于双通道测向系统的深度学习数据集生成方法
CN107202975A (zh) 一种二维矢量阵阵元姿态误差校正方法
CN113923590B (zh) 一种锚节点位置不确定情况下的toa定位方法
CN110673088A (zh) 混合视距和非视距环境中基于到达时间的目标定位方法
CN110221245A (zh) 联合估计目标位置和非视距误差的鲁棒tdoa定位方法
CN112835020B (zh) 面向非视距参数估计的刚体定位方法
CN112954637B (zh) 一种锚节点位置不确定情况下的目标定位方法
Jia-Jia et al. An effective frequency-spatial filter method to restrain the interferences for active sensors gain and phase errors calibration
CN112579972B (zh) 方向性电磁耦合效应下空域信息联合估计方法
CN107484119B (zh) 一种用于移动通信系统的终端跟踪定位方法
Wei et al. Angle–of–Arrival (AoA) Factorization in Multipath Channels

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