CN112180320B - 一种无人机无源定位系统及方法 - Google Patents

一种无人机无源定位系统及方法 Download PDF

Info

Publication number
CN112180320B
CN112180320B CN202010845680.2A CN202010845680A CN112180320B CN 112180320 B CN112180320 B CN 112180320B CN 202010845680 A CN202010845680 A CN 202010845680A CN 112180320 B CN112180320 B CN 112180320B
Authority
CN
China
Prior art keywords
time difference
order
difference extraction
value
signal
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
Application number
CN202010845680.2A
Other languages
English (en)
Other versions
CN112180320A (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.)
Chengdu University
Original Assignee
Chengdu 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 Chengdu University filed Critical Chengdu University
Priority to CN202010845680.2A priority Critical patent/CN112180320B/zh
Publication of CN112180320A publication Critical patent/CN112180320A/zh
Application granted granted Critical
Publication of CN112180320B publication Critical patent/CN112180320B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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
    • G01S5/00Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/02Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using radio waves

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明提供了一种无人机无源定位系统及方法,属于定位技术领域,本发明根据不同的实际时差值,选择与之对应的最优的傅里叶变换阶数,并进行卡尔曼滤波,并采用非递归算法Chan利用时差提取部分所提供的三个时差值获得初始坐标,然后,将无人机的初始估计坐标和时差值传送到递归算法Taylor部分进行迭代,用以校准误差,最终解算出高精度位置信息。解决了在利用FPGA实现傅里叶变换功能时,变换阶数的不合理选择会造成大量的资源浪费,同时,在两路信号的时差值与傅里叶变换阶数不耦合的情况下,将导致时差值提取出错的问题。

Description

一种无人机无源定位系统及方法
技术领域
本发明属于定位技术领域,尤其涉及一种无人机无源定位系统。
背景技术
目前国内外的无源定位技术,主要分为依靠非合作外辐射源的无源探测定位技术和依靠目标自身辐射源的无源定位技术,并且都得到了一定的发展和应用。但是,仍没有精度比较高的无源定位系统,同时,在功能方面,往往以追求高性能为主,较少考虑成本、数量等因素,因此,设计硬件资源可优化、时差提取精度较高的无源定位系统是非常有必要的。
无源定位技术早期主要应用在军事领域,在西方发达国家得到了重点研究应用,并取得了实战应用,例如:捷克ERA公司开发的VERA-E系统、俄罗斯无源定位系统卡尔秋塔、以色列EL-L8300G无源雷达系统等。近年来,随着“非合作型”无人机的广泛应用,对公共安全、保密场所的威胁日益增加。“非合作型”无人机一般为“低、小、慢”的目标,所以,通常采用基于TDOA(到达时间差) 体制的无人机无源定位技术进行定位,基于TDOA体制的时差提取算法主要包括基本互相关算法、广义互相关算法、自相关与互相关结合的算法和基于LMS自适应与广义互相关的结合算法等。
其中,基本互相关算法是一种最基本的算法,其算法原理简单,计算复杂程度低,计算量较小;但是,基本互相关算法要求信号与噪声、噪声与噪声之间不相关,而在实际环境中,可能会出现信号与噪声、噪声与噪声相关的情况,不利于工程实现;广义互相关算法(GCC)是由Knapp C和Carter G提出的,通过对信号进行加权滤波处理来提高信噪比,其主要的加权滤波器有ROTH处理器、平滑相干变换SCOT以及相位变换PHAT,但是,广义互相关算法,需要合理地选择加权函数,才能取得更好的滤波效果和更高的时差提取精度;LMS 算法是自适应地调节滤波器的权系数,以两路信号的均方误差为最小判据,进行时差提取,但是,LMS算法牺牲了计算速度来降低对信号和噪声统计先验知识的要求。更重要的是,以上算法都未考虑,在利用FPGA实现傅里叶变换功能时,变换阶数的不合理选择会造成大量的资源浪费,同时,在两路信号的时差值与傅里叶变换阶数不耦合的情况下,将导致时差值提取出错的问题。
发明内容
针对现有技术中的上述不足,本发明提供的一种无人机无源定位系统,解决了在利用FPGA实现傅里叶变换功能时,变换阶数的不合理选择会造成大量的资源浪费,同时,在两路信号的时差值与傅里叶变换阶数不耦合的情况下,将导致时差值提取出错的问题。
为了达到以上目的,本发明采用的技术方案为:
本方案提供一种无人机无源定位系统,包括依次连接的信号接收模块、时差提取模块以及定位解算模块;
所述信号接收模块,用于利用射频天线阵列接收无人机辐射出的无人机图传信号,并将所述无人机图传信号转换成数字信号;
所述时差提取模块,用于利用卡尔曼-最优阶互相关算法对所述数字信号进行时差提取处理,得到最优时差提取值;
所述定位解算模块,用于根据所述最优时差提取值,利用非递归算法计算得到无人机的初始估计坐标,并利用递归算法对无人机的初始估计坐标和最优时差提取值进行误差计算得到定位解算结果,完成对无人机的定位。
进一步地,所述时差提取模块包括快速傅里叶变换单元、互功率谱单元、快速傅里叶逆变换单元以及时延估计单元;
所述快速傅里叶变换单元,用于对数字信号进行归一化处理,得到频域信号,并对所述频域信号进行自相关处理,生成自功率频谱信号;
所述互功率谱单元,用于根据所述频域信号计算得到互功率频谱函数,以及根据自功率频谱信号计算得到二次相关功率频谱函数;
所述快速傅里叶逆变换单元,用于根据所述互功率频谱函数和二次相关功率频谱函数,计算得到广义互相关函数;
所述时延估计单元,用于取所述广义互相关函数的最大值时对应的自变量为快速傅里叶变换某一阶数内输出的时差提取值,并根据所述快速傅里叶变换某一阶数内输出的时差提取值计算得到最优时差提取值。
基于上述系统,本发明还提供了一种无人机无源定位方法,包括以下步骤:
S1、利用射频天线阵列接收无人机辐射出的无人机图传信号,并将所述无人机图传信号转换成数字信号;
S2、利用卡尔曼-最优阶互相关算法对所述数字信号进行时差提取处理,得到最优时差提取值;
S3、根据所述最优时差提取值,利用非递归算法计算得到无人机的初始估计坐标,并利用递归算法对无人机的初始估计坐标和最优时差提取值进行误差计算;
S4、判断所述误差计算结果是否小于预设的阈值,若是,则根据误差值输出定位解算结果,完成对无人机的定位,否则,返回步骤S3。
进一步地,所述步骤S2包括以下步骤:
S201、对所述数字信号进行幅值归一化处理,去除数字信号中的幅度信息;
S202、对归一化后的数字信号进行快速傅里叶变换,得到频域信号,并对所述频域信号进行自相关处理,生成自功率频谱信号;
S203、根据所述频域信号计算得到互功率频谱函数,以及根据自功率频谱信号计算得到二次相关功率频谱函数;
S204、根据所述互功率频谱函数和二次相关功率频谱函数,计算得到广义互相关函数;
S205、取所述广义互相关函数的最大值时对应的自变量为快速傅里叶变换某一阶数内输出的时差提取值,并根据所述快速傅里叶变换某一阶数内输出的时差提取值计算得到当前时刻对应阶数的时差提取值;
S206、根据所述当前时刻对应阶数的时差提取值,分别计算得到当前阶数与上一阶数的时差提取值误差以及当前阶数与下一阶数的时差提取值误差;
S207、当当前阶数与上一阶数的时差提取值误差不等于零,且当当前阶数与下一阶数的时差提取值误差等于零时,计算得到快速傅里叶变换最优阶数的时差估计值;
S208、根据所述快速傅里叶变换最优阶数的时差估计值,利用卡尔曼滤波计算得到最优时差提取值。
再进一步地,所述步骤S201中归一化处理的表达式如下:
Ci(t)=Xi(t)/max(abs(Xi(t)))
其中,Ci(t)表示归一化处理后的信号,Xi(t)表示接收的第i路信号,abs(·) 表示取绝对值,max(·)表示取数字信号中的最大值,i表示信号路数,且i=1,2,3,4。
再进一步地,所述步骤S202中自功率频谱信号的表达式如下:
Rii(ω)=Ri(ω)*conj(Ri(ω))
Ri(ω)=FFT(Ci(t),2N)
其中,Rii(ω)表示自功率频谱信号,Ri(ω)表示频域信号,conj(·)表示取共轭复数,Ci(t)表示归一化处理后的信号,FFT(·)表示快速傅里叶变换,N表示傅里叶变换阶数,i表示信号路数,且i=1,2,3,4。
再进一步地,所述步骤S204中广义互相关函数的表达式如下:
Yi,i+1(t)=IFFT(Gi,i+1(ω)*conj(Ki(ω)))
Gi,i+1(ω)=Ri(ω)*conj(Ri+1(ω))
Ki(ω)=Rii(ω)*conj(Ri+1,i+1(ω))
其中,Yi,i+1(t)表示广义互相关函数,IFFT(·)表示快速傅里叶逆变换,Gi,i+1(ω)表示互功率频谱函数,Ki(ω)表示二次相关功率频谱函数,conj(·)表示取共轭复数,i表示信号路数,且i=1,2,3,4,Ri(ω)表示频域信号,Ri+1(ω)表示第i+1路频谱信号,Rii(ω)表示自功率频谱信号,Ri+1,i+1(ω)表示第i+1路信号的自功率频谱信号,Yi,i+1(t)表示广义互相关函数。
再进一步地,所述步骤S205中当前时刻对应阶数的时差提取值的表达式如下:
Figure BDA0002642971500000051
Figure BDA0002642971500000052
Figure BDA0002642971500000053
其中,δn表示当前时刻对应阶数的时差提取值,
Figure BDA0002642971500000054
表示广义互相关函数取得最大值时对应的自变量,
Figure BDA0002642971500000055
表示修正参数,e表示指数,n表示最佳傅里叶变换阶数,mt-1表示由上一时刻的时差预测值校正后的最优时差提取值,mt-2表示t-2时刻的时差预测值校正后的最优时差提取值。
再进一步地,所述步骤S206中当前阶数与上一阶数的时差提取值误差的表达式如下:
θ1=δnn-1
其中,θ1表示当前阶数与上一阶数的时差提取值误差,δn表示当前时刻对应阶数的时差提取值,δn-1表示上一阶数的时差提取值;
所述当前阶数与下一阶数的时差提取值误差的表达式如下:
θ2=δnn+1
其中,θ2表示当前阶数与下一阶数的时差提取值误差,δn+1表示下一阶数的时差提取值。
再进一步地,所述步骤S208中最优时差提取值的表达式如下:
mt=mnew.t+Kt*(δn-mnew.t)
Figure BDA0002642971500000061
Pnew.t=Pt-1+Q
mnew.t=mt-1
其中,mt表示最优时差提取值,mnew.t表示当前时刻的时差预测值,Kt表示卡尔曼的计算增益,δn表示最优阶时差提取值,Pnew.t表示当前时刻最优时差估计值的方差,R表示两次测量时差值间的方差,Pt-1表示上一时刻的最优时差估计值的方差,Q表示连续两个时刻最优时差估计值的方差,mt-1表示由上一时刻的时差预测值校正后的最优时差提取值。
本发明的有益效果:
(1)本发明根据不同的实际时差值,选择与之对应的最优的傅里叶变换阶数,并进行卡尔曼滤波,既减少了对FPGA资源的占用,也避免了因傅里叶变换阶数导致的时差提取出错,同时,可以提高低信噪比条件下的抗噪能力时差提取精度。
(2)本发明的卡尔曼-最优阶互相关算法是将自适应获取最优傅里叶变换阶数与卡尔曼滤波两者进行结合的算法。卡尔曼-最优阶互相关算法的优点是通过对傅里叶变换阶数的扫描,可以使系统的硬件资源得以合理分配,同时,可以解决因傅里叶变换阶数导致的时差提取出错的问题。此外,卡尔曼滤波加快了最优时差提取的收敛速度,提升了时差提取的精度。
(3)本发明在定位解算部分采用由Chan和Taylor两部分定位算法。Chan 算法利用时差提取部分所提供的三个时差值获得初始坐标,然后,将初始估计坐标和时差值传送到Taylor算法部分进行迭代,用以校准误差,最终解算出高精度位置信息。
(4)本发明可解决当前国内外无人机无源定位系统实现成本高、硬件资源分配不合理、由阶数不耦合导致时差提取出差、时差提取精度不高和运行处理速度慢等问题;较低的成本,可使无源定位系统得以广泛的应用;资源合理分配,大大提高了系统的运行和处理速度;阶数的耦合,保证时差的准确提取;卡尔曼滤波,提高了低信噪比环境下的时差提取精度。
(5)本发明采用高效的、性能优异的时差提取技术,在机场、核电站、化工厂、政府大楼、军事基地、监狱等重要安防要地,对“非合作型”无人机的入侵防控,有着极为广泛的应用前景和实用空间。与此同时,本发明在警用作业、电力勘测、抗震救灾、野外临时通信保障和部分军事活动等领域还有着极为广泛的应用,市场空间巨大,经济效益可观。
附图说明
图1为本发明的系统结构示意图。
图2为本实施例中GUI人机交互界面设计图。
图3为本实施例中Y型基站布局图。
图4为本发明的方法流程图。
图5为本实施例中未加噪声的信号波形图。
图6为本实施例中阶数变换估计结果示意图。
图7为本实施例中参数修正结果示意图。
图8为本实施例中最优阶数时差值的变化结构示意图。
图9为本实施例中卡尔曼滤波结果示意图。
图10为本实施例中误差分析示意图。
图11为本实施例中五种相关算法曲线拟合性能比较图。
图12为本实施例中不同处理平台的数据处理速度对比图。
具体实施方式
下面对本发明的具体实施方式进行描述,以便于本技术领域的技术人员理解本发明,但应该清楚,本发明不限于具体实施方式的范围,对本技术领域的普通技术人员来讲,只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内,这些变化是显而易见的,一切利用本发明构思的发明创造均在保护之列。
实施例1
如图1所示,本发明提供了一种无人机无源定位系统,包括依次连接的信号接收模块、时差提取模块以及定位解算模块;
所述信号接收模块,用于利用射频天线阵列接收无人机辐射出的无人机图传信号,并将所述无人机图传信号转换成数字信号;
所述时差提取模块,用于利用卡尔曼-最优阶互相关算法对所述数字信号进行时差提取处理,得到最优时差提取值;
所述定位解算模块,用于根据所述最优时差提取值,利用非递归算法计算得到无人机的初始估计坐标,并利用递归算法对无人机的初始估计坐标和最优时差提取值进行误差计算得到定位解算结果,完成对无人机的定位。所述时差提取模块包括快速傅里叶变换单元、互功率谱单元、快速傅里叶逆变换单元以及时延估计单元;所述快速傅里叶变换单元,用于对数字信号进行归一化处理,得到频域信号,并对所述频域信号进行自相关处理,生成自功率频谱信号;所述互功率谱单元,用于根据所述频域信号计算得到互功率频谱函数,以及根据自功率频谱信号计算得到二次相关功率频谱函数;所述快速傅里叶逆变换单元,用于根据所述互功率频谱函数和二次相关功率频谱函数,计算得到广义互相关函数;所述时延估计单元,用于取所述广义互相关函数的最大值时对应的自变量为快速傅里叶变换某一阶数内输出的时差提取值,并根据所述快速傅里叶变换某一阶数内输出的时差提取值计算得到最优时差提取值。
本实施例中,如图2所示,本发明还包括GUI人机交互模块:GUI人机交互模块由MATLAB编写,在PC端对无人机的实时数据进行显示,内容包括四路基站的实时信号、载入的虚拟地图、真实的时差值、最终坐标相对位置以及最终坐标值的显示等功能。有利于测试者掌控和把握此时无人机的实时状态,同时,有利于整套系统的调试,也有利于对无人机的反制采取进一步的措施。
本实施例中,信号接收模块中的射频天线阵列组成的信号接收系统使用四单元天线阵列,单元宽带增益高,对于S、C波段的无人机信号,其增益高达 20dBi,对无人机常用频段均可支持,同时,使用面向3G和4G应用的高性能、高度集成的AD9361射频接收芯片,具备宽动态、低NF、高抗干扰、多通道一致性调节的能力。
本实施例中,时差提取模块的硬件设计依靠FPGA实现,基于含有ARM Cortex-a9的Xilinx Zynq-7000全可编程片上系统的嵌入式处理器,实现了硬件加速,使定位系统具有较高的实时性与刷新率。本实施例中,系统的数据交互通过GPIO完成FPGA的PL与PS间的交互,较传统FPGA+ARM的方式,通信速度更快,信息传递结构更简单。本系统具有便携性,可适应外场的工作环境,具有良好的环境适应性,适用于车载或地面站等设备。
本实施例中,快速傅里叶变换FFT单元是对四路具有相关性的信号作快速傅里叶变化,在该单元中,采用流式处理I/O(STREAMING I/O)设计模式,用于保证数据的连续输入和输出,流式处理I/O(STREAMING I/O)设计模式相对于基—2突发I/O结构、基—4突发I/O结构,其数据吞吐量最大,处理的数据最多。对于流式处理I/O结构,在处理当前帧数据变换的同时,可以加载下一帧输入数据,并输出前一帧的变换结果数据,支持连续的输入数据并连续的输出计算结果,在数据的输入端,传入数据的位宽是16比特,采用定点输入,并且不进行缩放。对于输入信号的控制,采用AXI总线协议,当发送确认TVALID和发送准备TREADY信号同时拉高时,输入数据为有效数据。驱动快速傅里叶变换FFT 单元必经数据缓存发送使能信号CONFIG_VALID,用于通知快速傅里叶变换 FFT单元准备接收数据。发送确认TVALID和发送准备TREADY是一对握手信号,应利用数据缓存先发送TVALID信号,快速傅里叶变换FFT单元再发送 TREADY信号。若数据缓存在快速傅里叶变换FFT单元发送TREADY信号后发送TVALID信号,容易出现锁死的情况,导致输出结果出错。每一个快速傅里叶变换FFT单元只处理一路信号,进行FFT变换的点数满足2的N次方。在数据的输出端,FFT蝶形运算变换后,输出数据的形式是以二进制倒序的形式排列,因此,采用自然顺序输出。为保证数据的精度,选择无压缩的数据输出模式,为了使输出数据的流畅、无误差传输,互功率谱模块须告诉快速傅里叶变换FFT单元已准备好接收数据,此时,快速傅里叶变换FFT单元才输出数据。
本实施例中,互功率谱单元用于求取FFT变换后两路数据的互功率频谱函数。在作FFT变换后,每一路数据的位宽以29比特的大小输出,在此单元中,首先需要对其中的一路信号进行取共轭复数变换,然后,再对两路信号进行复数乘法。
本实施例中,在进行IFFT之前,对信号进行了加权滤波和截位。互功率频谱单元输出的信号的位宽为55比特,而IFFT单元输入的数据位宽最大为34比特,因此,需要对信号进行截位。在截位时,考虑到低位信号可能被噪声所污染,故截去部分低位信号,保留高位信号。首先,找到最大的数据的位宽,保留好符号位,截去多余的符号位与部分低位数据,最后,使输入的数据位宽位 34比特。在该单元中,采用流式处理I/O(STREAMING I/O)设计模式,对于输入信号的控制,采用AXI总线协议,点数为2的N次方。
本实施例中,时延估计单元是对IFFT单元输出的信号序列进行时延估计。先将信号序列内的数据全部转换为正数,然后,找到序列中的最大值,最终,同时输出序列中最大值及其对应的时间,即为时差输出值。
本实施例的工作原理:首先,载入地图并点击开始,系统开始运行;基站开始接收无人机的信号,并将信号显示到GUI界面,同时,可以通过GUI界面检测系统是否正常运行,天线是否正常接收信号,其中,如图3所示,四单元射频天线阵列采用Y型布站方式,接收空中无人机辐射出的无人机图传信号 OFDM。将接收到的无人机图传信号OFDM由AD9361射频接收芯片换为数字信号,经统计时分复用技术传送至基于尔曼-最优阶互相关算法的时差提取模块进行时差解算。提取出的时差由AXI总线传至ARM端进行定位解算,利用Chan 算法获得无人机的初始估计坐标,并利用Taylor算法对无人机的初始估计坐标和时差值传进行迭代,将输出结果与阈值比较,进行定位更正,用以校准误差,满足条件时,输出最佳解算坐标,解算结果通过UART发送至GUI人机交互界面,完成坐标值及相对位置显示。
本实施例中,本发明在信道数据传输方面,使用统计时分复用技术。统计时分复用使用STDM帧来传送复用的数据,STDM帧通过按需动态地分配时隙,因此,统计时分复用可以提高线路的利用率,减少片上资源的消耗,同时,能够保证信号在处理时的连续性。本系统采用高效的、性能优异的时差提取技术,在机场、核电站、化工厂、政府大楼、军事基地、监狱等重要安防要地,对“非合作型”无人机的入侵防控,有着极为广泛的应用前景和实用空间。与此同时,本系统在警用作业、电力勘测、抗震救灾、野外临时通信保障和部分军事活动等领域还有着极为广泛的应用,市场空间巨大,经济效益可观。
实施例2
如图4所示,本发明还提供了一种无人机无源定位方法,其实现方法如下:
S1、利用射频天线阵列接收无人机辐射出的无人机图传信号,并将所述无人机图传信号转换成数字信号;
本实施例中,首先,载入地图并点击开始,系统开始运行;基站开始接收无人机图传信号,并将无人机图传信号显示到GUI界面,同时,可以通过GUI 界面检测系统是否正常运行,天线是否正常接收信号;接收到的模拟信号转换为数字信号后,经统计时分复用技术传送至FPGA并利用卡尔曼-最优阶互相关算法,进行时差提取。
本实施例中,射频天线阵列为四单元射频天线阵列,且其以Y型的方式进行布站。
S2、利用卡尔曼-最优阶互相关算法对所述数字信号进行时差提取处理,得到最优时差提取值,其实现方法如下:
S201、对所述数字信号进行幅值归一化处理,去除数字信号中的幅度信息;
本实施例中,将接收到的信号由AD9361射频接收芯片与FPGA的板载接口FMC传入,改进的广义互相关算法将两路信号进行幅值归一化,得到的结果为:
Ci(t)=Xi(t)/max(abs(Xi(t)))
其中,Ci(t)表示归一化处理后的信号,Xi(t)表示接收的第i路信号,abs(·) 表示取绝对值,max(·)表示取数字信号中的最大值,i表示信号路数,且i=1,2,3,4。
本实施例中,归一化去除了信号的幅度信息,只保留了信号的相位特性,对于噪声和混响都有较好的抑制作用。
S202、对归一化后的数字信号进行快速傅里叶变换,得到频域信号,并对所述频域信号进行自相关处理,生成自功率频谱信号,其表达式如下:
Rii(ω)=Ri(ω)*conj(Ri(ω))
Ri(ω)=FFT(Ci(t),2N)
其中,Rii(ω)表示自功率频谱信号,Ri(ω)表示频域信号,conj(·)表示取共轭复数,Ci(t)表示归一化处理后的信号,FFT(·)表示快速傅里叶变换,N表示傅里叶变换阶数。
S203、根据所述频域信号计算得到互功率频谱函数,以及根据自功率频谱信号计算得到二次相关功率频谱函数;
S204、根据所述互功率频谱函数和二次相关功率频谱函数,计算得到广义互相关函数,其表达式如下:
Yi,i+1(t)=IFFT(Gi,i+1(ω)*conj(Ki(ω)))
Gi,i+1(ω)=Ri(ω)*conj(Ri+1(ω))
Ki(ω)=Rii(ω)*conj(Ri+1,i+1(ω))
其中,Yi,i+1(t)表示广义互相关函数,IFFT(·)表示快速傅里叶逆变换,Gi,i+1(ω)表示互功率频谱函数,Ki(ω)表示二次相关功率频谱函数,conj(·)表示取共轭复数,Ri(ω)表示频域信号,Ri+1(ω)表示第i+1路频谱信号,Rii(ω)表示自功率频谱信号,Ri+1,i+1(ω)表示第i+1路信号的自功率频谱信号。
S205、取所述广义互相关函数的最大值时对应的自变量为快速傅里叶变换某一阶数内输出的时差提取值,并根据所述快速傅里叶变换某一阶数内输出的时差提取值计算得到当前时刻对应阶数的时差提取值,其表达式如下:
Figure BDA0002642971500000141
Figure BDA0002642971500000142
Figure BDA0002642971500000143
其中,δn表示当前时刻对应阶数的时差提取值,
Figure BDA0002642971500000144
表示广义互相关函数取得最大值时对应的自变量,
Figure BDA0002642971500000145
表示修正参数,e表示指数,n表示最佳傅里叶变换阶数,mt-1表示由上一时刻的时差预测值校正后的最优时差提取值,mt-2表示 t-2时刻的时差预测值校正后的最优时差提取值。
S206、根据所述当前时刻对应阶数的时差提取值,分别计算得到当前阶数与上一阶数的时差提取值误差以及当前阶数与下一阶数的时差提取值误差;
S207、当当前阶数与上一阶数的时差提取值误差不等于零时,且当当前阶数与下一阶数的时差提取值误差等于零时,计算得到快速傅里叶变换最优阶数的时差估计值;
本实施例中,当进行时差值提取的快速傅里叶变换时,提取出的时差值随快速傅里叶变换的阶数而改变,但是,当达到一定的阶数后,提取出的时差保持正确且不再改变,称这一阶数为最优阶。
令:
θ1=δnn-1
θ2=δnn+1
当θ1≠0且θ2=0时,n为最佳傅里叶变换阶数,θ1表示当前阶数与上一阶数的时差提取值误差,θ2表示当前阶数与下一阶数的时差提取值误差。
Figure BDA0002642971500000151
其中,
Figure BDA0002642971500000152
表示快速傅里叶变换最优阶数的时差估计值。
S207、根据当前阶数与上一阶数的时差提取值误差以及当前阶数与下一阶数的时差提取值误差,计算得到快速傅里叶变换最优阶数的时差估计值;
S208、根据所述快速傅里叶变换最优阶数的时差估计值,利用卡尔曼滤波计算得到最优时差提取值
本实施例中,将前四次最优阶数的时差估计值的平均值作为卡尔曼滤波的初始时差提取值m1
Figure BDA0002642971500000153
将上一时刻的最优估计值作为当前时刻的时差预测值:
mnew.t=mt-1
其中,mnew.t表示当前时刻的时差预测值,mt-1表示由上一时刻的时差预测值校正后的最优时差提取值。
将上一时刻最优时差估计值的方差加上过程方差作为当前时刻的方差预测值:
Pnew.t=Pt-1+Q
其中,Pnew.t表示当前时刻最优时差估计值的方差,Pi-1为上一时刻最优时差估计值的方差,Q表示连续两个时刻最优时差估计值的方差。
根据上一时刻最优时差估计值的方差计算得到卡尔曼的计算增益:
Figure BDA0002642971500000154
其中,Kt表示卡尔曼的计算增益,R表示两次测量时差值间的方差。
结合当前时刻时差的提取值,对上一时刻的预测进行校正,并输出最优时差提取值:
mt=mnew.t+Kt*(δn-mnew.t)
其中,mt表示最优时差提取值,mnew.t表示当前时刻的时差预测值,Kt表示卡尔曼的计算增益,δn表示最优阶时差提取值。
对最优时差估计值的方差进行更新:
Pi=(1-Ki)*Pnewi
其中,Pi表示当前时刻最优时差估计值的方差。
本实施例中,以卡尔曼—最优阶互相关算法为核心的时差提取算法的优点是通过对傅里叶变换阶数的扫描,可以使系统的硬件资源得以合理分配,同时,可以解决因傅里叶变换阶数导致的时差提取出错的问题。此外,卡尔曼滤波加快了最优时差提取的收敛速度,提升了时差提取的精度,对于无人机图传信号 (OFDM信号)展现出了更强的资源分配能力和抗噪能力。
S3、根据所述最优时差提取值,利用非递归算法Chan计算得到无人机的初始估计坐标,并利用递归算法Taylor对所述无人机的初始估计坐标和最优时差提取值进行误差计算;
S4、判断所述误差计算结果是否小于预设的阈值,若是,则根据误差值输出定位解算结果,完成定位,否则,返回步骤S3。
本实施例中,本申请提出Chan-Taylor联合算法,将Chan算法解算出的目标坐标作为无人机的初始估计坐标值赋给Taylor算法进行迭代运算,即使获取的时差值存在一定误差,令无人机的初始估计坐标的精度不高,但可以通过迭代提高定位坐标的精度。
为了进一步对本发明进行说明,现以以下实验数据进行说明。
通过仿真对比分析相同条件下尔曼-最优阶互相关算法和基本互相关算法、广义加权互相关算法的时差提取精度,验证其更优的性能。仿真和试验结果表明,相对于传统算法,卡尔曼-最优阶互相关算法在较低信噪比的情况下,具有更强的抗噪能力和更高的时差提取精度。仿真分析中使用的广义互相关加权函数如表1所示。
表1
Figure BDA0002642971500000171
表1中,ξ(f)表示为广义互相关频域的加权函数,Φ11(f)和Φ22(f)分别表示信号X1(t)和X2(t)的自相关函数,当ξ(f)=1时,即为基于互相关算法。
本实施例中,实验利用MATLAB进行仿真,仿真条件为:采集的某无人机图传信号,采样频率Fs=200MHZ,采样周期Ts=5ns。在仿真中,两段信号的时差为122Ts,信号波形如图5所示。在仿真中,选取时差值为122个点的两路信号,傅里叶变换的阶数以2的1次方为公比,从2的1次方变换到2的 20次方,实验结果如图6所示。
本实施例中,如图7所示,在时差值为122的条件下,当傅里叶变换的阶数为2的9次方以后,时差提取值保持稳定的输出结果,故9阶为最优傅里叶变换阶数。在此基础上,连续改变两路信号的时差值时,最优阶数的变换如图8 所示,在图8中,最优阶数随着时差值的增加而增大,说明两路信号的时差值与傅里叶变换阶数存在关系。将两值的关系统一起来,得到一个与时差变化相关的量δn,δn满足任一时差值与最优阶数的关系。将最优阶得到的时差值进行卡尔曼滤波,结果如图9所示,在图9中,在50μs内,卡尔曼-最优阶互相关算法具有较好的收敛值,可以更准确地估计出实际的时差值,并提高算法在不同信噪比条件下的时差提取精度。误差分析如图10所示,在图9、图10中,在 50μs内,卡尔曼-最优阶互相关算法具有较好的收敛值,可以更准确地估计出实际的时差值,并提高算法在不同信噪比条件下的时差提取精度。
本实施例中,对算法进行曲线拟合,比较五种算法的抗噪能力、时差提取精度和收敛速度,仿真结果如图11所示,通过图11对比分析可以发现,随着信噪比的降低,相对于其它四种算法,卡尔曼-最优阶互相关算法,无论是在抗躁能力和时差提取精度上,还是在收敛速度上,均有明显的优势,其中,ROTH 加权函数的误差最大,说明它对高斯白噪声滤波效果并不理想。可见,合适的加权函数是提高时差提取精度的重要因素。为测试系统的性能,在相同数据吞吐量条件下,对比不同处理平台的数据处理速度,不同平台的对比结果如图12 所示,由图12的实验结果,可以分析在不同平台的处理速度,其中,MATLAB 的处理时间为0.254s,C#的处理时间为0.0408362s,而FPGA的处理时间为579us。因此,本发明在实际工程应用上具有明显的优势。

Claims (3)

1.一种无人机无源定位系统,其特征在于,包括依次连接的信号接收模块、时差提取模块以及定位解算模块;
所述信号接收模块,用于利用射频天线阵列接收无人机辐射出的无人机图传信号,并将所述无人机图传信号转换成数字信号;
所述时差提取模块,用于利用卡尔曼-最优阶互相关算法对所述数字信号进行时差提取处理,得到最优时差提取值,其具体为:
A1、对所述数字信号进行幅值归一化处理,去除数字信号中的幅度信息;
A2、对归一化后的数字信号进行快速傅里叶变换,得到频域信号,并对所述频域信号进行自相关处理,生成自功率频谱信号;
A3、根据所述频域信号计算得到互功率频谱函数,以及根据自功率频谱信号计算得到二次相关功率频谱函数;
A4、根据所述互功率频谱函数和二次相关功率频谱函数,计算得到广义互相关函数;
A5、取所述广义互相关函数的最大值时对应的自变量为快速傅里叶变换某一阶数内输出的时差提取值,并根据所述快速傅里叶变换某一阶数内输出的时差提取值进行最优阶参数修正计算得到当前时刻对应阶数的时差提取值;
A6、根据所述当前时刻对应阶数的时差提取值,分别计算得到当前阶数与上一阶数的时差提取值误差以及当前阶数与下一阶数的时差提取值误差;
A7、当当前阶数与上一阶数的时差提取值误差不等于零,且当前阶数与下一阶数的时差提取值误差等于零时,计算得到快速傅里叶变换最优阶数的时差估计值;
A8、根据所述快速傅里叶变换最优阶数的时差估计值,利用卡尔曼滤波计算得到最优时差提取值;
所述步骤A1中归一化处理的表达式如下:
Ci(t)=Xi(t)/max(abs(Xi(t)))
其中,Ci(t)表示归一化处理后的信号,Xi(t)表示接收的第i路信号,abs(·)表示取绝对值,max(·)表示取数字信号中的最大值,i表示信号路数,且i=1,2,3,4;
所述步骤A2中自功率频谱信号的表达式如下:
Rii(ω)=Ri(ω)*conj(Ri(ω))
Ri(ω)=FFT(Ci(t),2N)
其中,Rii(ω)表示自功率频谱信号,Ri(ω)表示频域信号,conj(·)表示取共轭复数,Ci(t)表示归一化处理后的信号,FFT(·)表示快速傅里叶变换,N表示快速傅里叶变换阶数,i表示信号路数,且i=1,2,3,4;
所述步骤A4中广义互相关函数的表达式如下:
Yi,i+1(t)=IFFT(Gi,i+1(ω)*conj(Ki(ω)))
Gi,i+1(ω)=Ri(ω)*conj(Ri+1(ω))
Ki(ω)=Rii(ω)*conj(Ri+1,i+1(ω))
其中,Yi,i+1(t)表示广义互相关函数,IFFT(·)表示快速傅里叶逆变换,Gi,i+1(ω)表示互功率频谱函数,Ki(ω)表示二次相关功率频谱函数,conj(·)表示取共轭复数,i表示信号路数,且i=1,2,3,4,Ri(ω)表示频域信号,Ri+1(ω)表示第i+1路频谱信号,Rii(ω)表示自功率频谱信号,Ri+1,i+1(ω)表示第i+1路信号的自功率频谱信号;
所述步骤A5中当前时刻对应阶数的时差提取值的表达式如下:
Figure FDA0002968811490000031
Figure FDA0002968811490000032
Figure FDA0002968811490000033
其中,δn表示当前时刻对应阶数的时差提取值,
Figure FDA0002968811490000034
表示广义互相关函数取得最大值时对应的自变量,
Figure FDA0002968811490000035
表示修正参数,e表示指数,n表示最佳傅里叶变换阶数,mt-1表示由上一时刻的时差预测值校正后的最优时差提取值,mt-2表示t-1上一时刻的时差预测值校正后的最优时差提取值,Yi,i+1(t)表示广义互相关函数;
所述步骤A6中当前阶数与上一阶数的时差提取值误差的表达式如下:
θ1=δnn-1
其中,θ1表示当前阶数与上一阶数的时差提取值误差,δn表示当前时刻对应阶数的时差提取值,δn-1表示上一阶数的时差提取值;
所述当前阶数与下一阶数的时差提取值误差的表达式如下:
θ2=δnn+1
其中,θ2表示当前阶数与下一阶数的时差提取值误差,δn+1表示下一阶数的时差提取值;
所述步骤A8中最优时差提取值的表达式如下:
mt=mnew.t+Kt*(δn-mnew.t)
Figure FDA0002968811490000036
Pnew.t=Pt-1+Q
mnew.t=mt-1
其中,mt表示最优时差提取值,mnew.t表示当前时刻的时差预测值,Kt表示卡尔曼的计算增益,δn表示最优阶时差提取值,Pnew.t表示当前时刻最优时差估计值的方差,R表示两次测量时差值间的方差,Pt-1表示上一时刻的最优时差估计值的方差,Q表示连续两个时刻最优时差估计值的方差,mt-1表示由上一时刻的时差预测值校正后的最优时差提取值;
所述定位解算模块,用于根据所述最优时差提取值,利用非递归算法计算得到无人机的初始估计坐标,并利用递归算法对无人机的初始估计坐标和最优时差提取值进行误差计算得到定位解算结果,完成对无人机的定位。
2.根据权利要求1所述的无人机无源定位系统,其特征在于,所述时差提取模块包括快速傅里叶变换单元、互功率谱单元、快速傅里叶逆变换单元以及时延估计单元;
所述快速傅里叶变换单元,用于对数字信号进行归一化处理,得到频域信号,并对所述频域信号进行自相关处理,生成自功率频谱信号;
所述互功率谱单元,用于根据所述频域信号计算得到互功率频谱函数,以及根据自功率频谱信号计算得到二次相关功率频谱函数;
所述快速傅里叶逆变换单元,用于根据所述互功率频谱函数和二次相关功率频谱函数,计算得到广义互相关函数;
所述时延估计单元,用于取所述广义互相关函数的最大值时对应的自变量为快速傅里叶变换某一阶数内输出的时差提取值,并根据所述快速傅里叶变换某一阶数内输出的时差提取值计算得到最优时差提取值。
3.一种无人机无源定位方法,其特征在于,包括以下步骤:
S1、利用射频天线阵列接收无人机辐射出的无人机图传信号,并将所述无人机图传信号转换成数字信号;
S2、利用卡尔曼-最优阶互相关算法对所述数字信号进行时差提取处理,得到最优时差提取值;
S3、根据所述最优时差提取值,利用非递归算法计算得到无人机的初始估计坐标,并利用递归算法对无人机的初始估计坐标和最优时差提取值进行误差计算;
S4、判断所述误差计算结果是否小于预设的阈值,若是,则根据误差值输出定位解算结果,完成对无人机的定位,否则,返回步骤S3;
所述步骤S2包括以下步骤:
S201、对所述数字信号进行幅值归一化处理,去除数字信号中的幅度信息;
S202、对归一化后的数字信号进行快速傅里叶变换,得到频域信号,并对所述频域信号进行自相关处理,生成自功率频谱信号;
S203、根据所述频域信号计算得到互功率频谱函数,以及根据自功率频谱信号计算得到二次相关功率频谱函数;
S204、根据所述互功率频谱函数和二次相关功率频谱函数,计算得到广义互相关函数;
S205、取所述广义互相关函数的最大值时对应的自变量为快速傅里叶变换某一阶数内输出的时差提取值,并根据所述快速傅里叶变换某一阶数内输出的时差提取值进行最优阶参数修正计算得到当前时刻对应阶数的时差提取值;
S206、根据所述当前时刻对应阶数的时差提取值,分别计算得到当前阶数与上一阶数的时差提取值误差以及当前阶数与下一阶数的时差提取值误差;
S207、当当前阶数与上一阶数的时差提取值误差不等于零,且当前阶数与下一阶数的时差提取值误差等于零时,计算得到快速傅里叶变换最优阶数的时差估计值;
S208、根据所述快速傅里叶变换最优阶数的时差估计值,利用卡尔曼滤波计算得到最优时差提取值;
所述步骤S201中归一化处理的表达式如下:
Ci(t)=Xi(t)/max(abs(Xi(t)))
其中,Ci(t)表示归一化处理后的信号,Xi(t)表示接收的第i路信号,abs(·)表示取绝对值,max(·)表示取数字信号中的最大值,i表示信号路数,且i=1,2,3,4;
所述步骤S202中自功率频谱信号的表达式如下:
Rii(ω)=Ri(ω)*conj(Ri(ω))
Ri(ω)=FFT(Ci(t),2N)
其中,Rii(ω)表示自功率频谱信号,Ri(ω)表示频域信号,conj(·)表示取共轭复数,Ci(t)表示归一化处理后的信号,FFT(·)表示快速傅里叶变换,N表示快速傅里叶变换阶数,i表示信号路数,且i=1,2,3,4;
所述步骤S204中广义互相关函数的表达式如下:
Yi,i+1(t)=IFFT(Gi,i+1(ω)*conj(Ki(ω)))
Gi,i+1(ω)=Ri(ω)*conj(Ri+1(ω))
Ki(ω)=Rii(ω)*conj(Ri+1,i+1(ω))
其中,Yi,i+1(t)表示广义互相关函数,IFFT(·)表示快速傅里叶逆变换,Gi,i+1(ω)表示互功率频谱函数,Ki(ω)表示二次相关功率频谱函数,conj(·)表示取共轭复数,i表示信号路数,且i=1,2,3,4,Ri(ω)表示频域信号,Ri+1(ω)表示第i+1路频谱信号,Rii(ω)表示自功率频谱信号,Ri+1,i+1(ω)表示第i+1路信号的自功率频谱信号;
所述步骤S205中当前时刻对应阶数的时差提取值的表达式如下:
Figure FDA0002968811490000071
Figure FDA0002968811490000072
Figure FDA0002968811490000073
其中,δn表示当前时刻对应阶数的时差提取值,
Figure FDA0002968811490000074
表示广义互相关函数取得最大值时对应的自变量,
Figure FDA0002968811490000075
表示修正参数,e表示指数,n表示最佳傅里叶变换阶数,mt-1表示由上一时刻的时差预测值校正后的最优时差提取值,mt-2表示t-1上一时刻的时差预测值校正后的最优时差提取值,Yi,i+1(t)表示广义互相关函数;
所述步骤S206中当前阶数与上一阶数的时差提取值误差的表达式如下:
θ1=δnn-1
其中,θ1表示当前阶数与上一阶数的时差提取值误差,δn表示当前时刻对应阶数的时差提取值,δn-1表示上一阶数的时差提取值;
所述当前阶数与下一阶数的时差提取值误差的表达式如下:
θ2=δnn+1
其中,θ2表示当前阶数与下一阶数的时差提取值误差,δn+1表示下一阶数的时差提取值;
所述步骤S208中最优时差提取值的表达式如下:
mt=mnew.t+Kt*(δn-mnew.t)
Figure FDA0002968811490000076
Pnew.t=Pt-1+Q
mnew.t=mt-1
其中,mt表示最优时差提取值,mnew.t表示当前时刻的时差预测值,Kt表示卡尔曼的计算增益,δn表示最优阶时差提取值,Pnew.t表示当前时刻最优时差估计值的方差,R表示两次测量时差值间的方差,Pt-1表示上一时刻的最优时差估计值的方差,Q表示连续两个时刻最优时差估计值的方差,mt-1表示由上一时刻的时差预测值校正后的最优时差提取值。
CN202010845680.2A 2020-08-20 2020-08-20 一种无人机无源定位系统及方法 Expired - Fee Related CN112180320B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010845680.2A CN112180320B (zh) 2020-08-20 2020-08-20 一种无人机无源定位系统及方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010845680.2A CN112180320B (zh) 2020-08-20 2020-08-20 一种无人机无源定位系统及方法

Publications (2)

Publication Number Publication Date
CN112180320A CN112180320A (zh) 2021-01-05
CN112180320B true CN112180320B (zh) 2021-04-27

Family

ID=73924166

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010845680.2A Expired - Fee Related CN112180320B (zh) 2020-08-20 2020-08-20 一种无人机无源定位系统及方法

Country Status (1)

Country Link
CN (1) CN112180320B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113358931B (zh) * 2021-05-14 2022-08-23 深圳华创电科技术有限公司 一种基于互功率谱的时差计算方法
CN113484823B (zh) * 2021-06-21 2024-03-29 南京航空航天大学 一种基于闭式补偿的高分辨率时延估计方法
CN113572540B (zh) * 2021-06-28 2023-04-18 中国电子科技集团公司第三十八研究所 一种基于相关域检测的无人机图传信号识别方法及系统
CN114047476B (zh) * 2022-01-13 2022-04-08 中国人民解放军空军预警学院 一种基于无人机集群的无源定位方法及系统

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6407703B1 (en) * 2000-08-07 2002-06-18 Lockheed Martin Corporation Multi-platform geolocation method and system
JP3528808B2 (ja) * 2001-04-03 2004-05-24 日本電気株式会社 相互相関関数計算方法および装置
CN107340497B (zh) * 2017-07-06 2019-01-08 中国人民解放军火箭军研究院 一种基于频域互相关的分布式时差测量方法
CN109143147A (zh) * 2018-10-23 2019-01-04 烟台东方威思顿电气有限公司 基于卡尔曼滤波算法的电能表检定装置功率测量方法
CN110515038B (zh) * 2019-08-09 2023-03-28 达洛科技(广州)有限公司 一种基于无人机-阵列的自适应无源定位装置及实现方法

Also Published As

Publication number Publication date
CN112180320A (zh) 2021-01-05

Similar Documents

Publication Publication Date Title
CN112180320B (zh) 一种无人机无源定位系统及方法
CN111010249A (zh) 一种角度时延域信道预测方法、预测系统及应用
CN106992800B (zh) 基于迭代自适应算法的电力线通信系统脉冲噪声抑制方法
CN105137454B (zh) 一种基于协方差矩阵特征分解的抗干扰算法的fpga实现方法及实现装置
Wang et al. The robust sparse Fourier transform (RSFT) and its application in radar signal processing
CN113542162B (zh) 基于块稀疏贝叶斯算法的上下链路通信感知一体化方法
CN106027445A (zh) 一种水声块结构稀疏特性的信道估计方法
CN112272068B (zh) 一种基于多任务压缩感知的多样化干扰估计和抑制方法
CN102158443A (zh) 一种抑制多分量lfm信号时频分布中交叉项的方法
WO2023000614A1 (zh) 无线定位参数估计方法、装置、系统、计算机设备及存储介质
CN113259283B (zh) 一种基于循环神经网络的单通道时频混叠信号盲分离方法
Yang et al. Underdetermined DOA estimation for moving array with reduced mutual coupling in unknown nonuniform noise environment
CN104914451B (zh) 一种块Toeplitz矩阵低复杂度求逆的空时抗干扰方法
CN109412984B (zh) 一种基于Aitken加速法的多天线场景下盲信噪比估算方法
CN112383492A (zh) 应用于短波ofdm双选天波信道估计的递归压缩感知方法及系统
CN106908760A (zh) 基于阵列自相关矩阵的单站无源定位方法
CN104363078B (zh) 基于鲁棒竞争聚类的欠定系统实正交空时分组码盲识别方法
Changgan et al. An improved forward/backward spatial smoothing root-MUSIC algorithm based on signal decorrelation
CN113671477B (zh) 一种基于图信号处理的雷达目标距离估计方法
CN115052246A (zh) 一种未知衰减系数下基于多频率代价函数融合的宽带信号直接定位方法
CN109932681B (zh) 一种基于空-时信息的降冗余嵌套阵列设置方法
CN109787676B (zh) 一种高动态下的零陷展宽方法
Liu et al. A real-time, pipelined incoherent dedispersion method and implementation in FPGA
Das et al. Performance analysis of TLS-Esprit and QR TLS-Esprit algorithm for Direction of Arrival estimation
CN106338743B (zh) 基于ccd解算的空频联合自适应调零算法

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210427

CF01 Termination of patent right due to non-payment of annual fee