CN113075461B - 基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法 - Google Patents

基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法 Download PDF

Info

Publication number
CN113075461B
CN113075461B CN202110194898.0A CN202110194898A CN113075461B CN 113075461 B CN113075461 B CN 113075461B CN 202110194898 A CN202110194898 A CN 202110194898A CN 113075461 B CN113075461 B CN 113075461B
Authority
CN
China
Prior art keywords
signal
lightning
antenna
time
window
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
CN202110194898.0A
Other languages
English (en)
Other versions
CN113075461A (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.)
Zhuhai Fudan Innovation Research Institute
Original Assignee
Zhuhai Fudan Innovation Research Institute
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 Zhuhai Fudan Innovation Research Institute filed Critical Zhuhai Fudan Innovation Research Institute
Priority to CN202110194898.0A priority Critical patent/CN113075461B/zh
Priority to PCT/CN2021/082299 priority patent/WO2022174486A1/zh
Publication of CN113075461A publication Critical patent/CN113075461A/zh
Application granted granted Critical
Publication of CN113075461B publication Critical patent/CN113075461B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R29/00—Arrangements for measuring or indicating electric quantities not covered by groups G01R19/00 - G01R27/00
    • G01R29/08—Measuring electromagnetic field characteristics
    • G01R29/0807—Measuring electromagnetic field characteristics characterised by the application
    • G01R29/0814—Field measurements related to measuring influence on or from apparatus, components or humans, e.g. in ESD, EMI, EMC, EMP testing, measuring radiation leakage; detecting presence of micro- or radiowave emitters; dosimetry; testing shielding; measurements related to lightning
    • G01R29/0842—Measurements related to lightning, e.g. measuring electric disturbances, warning systems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R29/00—Arrangements for measuring or indicating electric quantities not covered by groups G01R19/00 - G01R27/00
    • G01R29/08—Measuring electromagnetic field characteristics
    • G01R29/0864—Measuring electromagnetic field characteristics characterised by constructional or functional features
    • G01R29/0878—Sensors; antennas; probes; detectors
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R29/00—Arrangements for measuring or indicating electric quantities not covered by groups G01R19/00 - G01R27/00
    • G01R29/08—Measuring electromagnetic field characteristics
    • G01R29/0864—Measuring electromagnetic field characteristics characterised by constructional or functional features
    • G01R29/0892—Details related to signal analysis or treatment; presenting results, e.g. displays; measuring specific signal features other than field strength, e.g. polarisation, field modes, phase, envelope, maximum value
    • 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
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/02—Position-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
    • G01S5/06—Position of source determined by co-ordinating a plurality of position lines defined by path-difference measurements
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • 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
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明属于闪电定位技术领域，具体为基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法。本发明结合干涉仪的超短基线优势和多站到达时间差技术的三维定位能力，设计超短基线多站闪电观测布局，摒弃传统的“质心法”技术路线，通过对原始信号进行增采样，引入集合经验模态分解技术对增采样信号进行分析和优化，采用基于主窗口与辅助窗口相结合的互相关信号匹配技术，实现对甚高频辐射信号中脉冲信息的充分匹配和提取，利用到达时间差技术进行三维定位。本发明可以实现对闪电放电过程超高分辨率的三维定位，定位闪电放电事件的时间分辨率可以达到纳秒量级，空间误差的理论分辨率提升到分米量级。

Description

基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法
技术领域
本发明属于闪电定位技术领域，具体涉及一种基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法。
背景技术
闪电发生时放射的不同频段的电磁信号均可用于闪电定位。根据探测系统的差异和定位方法的不同，闪电定位方法主要分为三类：磁定向法(Magnetic DirectionFinding,MDF)、到达时间差法(Time of Arrival,TOA)和干涉法(Interferometry)。
到达时间差的定位技术是通过测量特定频段闪电辐射的电磁信号到达不同辐射接收天线的时间差来计算辐射源发生的位置信息，根据接收天线布设距离可大致分为三类：(1)超短基线(几十到几百米)、(2)短基线(几公里至几十公里)、(3)长基线(几十至几百公里)。短基线探测系统接收闪电VHF频段(30-300MHz)信号；长基线系统接收闪电VLF和LF频段(3-300KHz)信号。目前国际上相关研究领域应用最广、影响最大的运行在VHF频段的TOA辐射定位系统是美国(New Mexico Institute of Mining and Technology)发展的闪电VHF辐射源定位系统(Lightning Mapping Array,LMA)(Thomasand Ronald,2004),通常由10个左右工作在60-66MHz的辐射接收天线组网观测，闪电辐射源定位的时间分辨率在几十微秒量级，空间定位误差为几十米。长基线定位系统与LMA的定位原理相似，特点是站网覆盖范围更广，但是时空分辨率通常都不高(Shao et al.,2006)。
随着各方面技术水平特别是电子和计算机技术的快速发展，促进了基于干涉法的闪电探测系统的发展。基于宽带甚高频信号探测的宽带干涉仪技术和闪电定位方法被Shaoet al.(1996)最早引入到闪电观测和研究领域。目前常见的宽带干涉仪的天线布局继承了窄带干涉仪系统基线间保持垂直的正交基线结构，即由三个VHF辐射接收天线构成两条正交的基线，基线长度一般几十米到上百米，天线的工作频带较宽，通常在几十MHz到几百MHz之间。通过探测闪电射频信号到达天线的时间差来实现对闪电放电过程的二维定位。这类系统由于基线非常短，通常采用单个高性能采集卡对三个天线数据同步采集，采样率高达几百兆，采样点的时间分辨率高达纳秒量级，定位得到的辐射源的时间分辨率可以达到亚微秒量级(Stocket al.,2014)。
但无论是目前主流的到达时间差法(TOA)还是干涉仪闪电定位系统，都有各自的优缺点：干涉仪的基线很短，不同天线间的信号非常容易被匹配求取时间差，因而定位结果的时间分辨率非常高，现有技术水平下可以达到亚微秒量级，不足之处在于干涉仪仅能获取闪电放电事件的二维信息(方位角和仰角)，并且由于采用平面波近似，其定位理论本身存在非常显著的系统误差；基于TOA的LMA类似的闪电定位系统虽然能够获得闪电放电过程的三维信息，但是由于基线长度相对较长，不同天线上的信号差异会很大，因而信号匹配的难度较大，准确率较低，这类系统目前能够达到的对闪电放电事件的时间分辨率最高为几十微秒量级，空间误差为几十米；此外，目前主流的干涉仪和三维闪电定位系统均普遍采用“质心法(centroid method)”来匹配信号求取时间差，这个技术实际上限制了VHF观测信号潜力的充分挖掘。
随着闪电物理研究的不断深入，对闪电定位系统获取的定位信息的时空分辨率的要求不断提升。基于现实的需求以及对现有多种闪电定位系统的优缺点的分析，本发明提出一种基于宽带甚高频辐射信号的超短基线闪电三维定位方法和系统，即结合干涉仪的超短基线优势和多站TOA技术的三维定位能力，提出一种超短基线多站观测布局，摒弃传统的“质心法”技术路线，通过对原始信号进行增采样，引入集合经验模态分解技术对增采样信号进行分析和优化，采用基于主窗口与辅助窗口相结合的互相关信号匹配技术，实现VHF辐射信号中脉冲信息的充分匹配和提取，利用TOA技术进行三维定位。该技术可以实现对闪电放电过程超高分辨率的三维定位，定位事件的时间分辨率可以达到纳秒量级，空间误差的理论分辨率提升到分米量级。
参考文献
Shao,X.M.,Stanley,M.,Regan,A.,Harlin,J.,Pongratz,M.,&Stock,M.(2006).Total lightning observations with the new and improved Los Alamos sfericarray(LASA).Journal of Atmospheric&Oceanic Technology,23(23),1273,doi:10.1175/JTECH1908.1.
Stock,M.G.,Akita,M.,Krehbiel,P.R.,Rison,W.,Edens,H.E.,Kawazaki,Z.,&Stanley,M.A.(2014).Continuous broadband digital interferometry of lightningusing a generalized cross-correlation algorithm.Journal of GeophysicalResearch:Atmospheres,119,3134–3165.https://doi.org/10.1002/2013JD020217
Stock,M.G.,P.R.Krehbiel,J.Lapierre,T.Wu,M.A.Stanley and H.E.Edens(2017),Fast positive breakdown in lightning,Journal of Geophysical Research:Atmospheres,122,15,8135-8152.
Thomas,&Ronald,J..(2004).Accuracy of the lightning mappingarray.Journal of Geophysical Research,109(D14),D14207.。
发明内容
本发明的目的在于提供一种时间和空间分辨率高的基于宽带甚高频辐射信号的超短基线闪电三维定位方法。
本发明提供的基于宽带甚高频辐射信号的超短基线闪电三维定位方法，结合干涉仪的超短基线定位优势和多站到达时间差技术的三维定位能力，采用超短基线多站观测布局；摒弃传统的“质心法”技术路线，通过对原始信号进行增采样；引入集合经验模态分解技术对增采样信号进行分析和优化；采用基于主窗口与辅助窗口相结合的互相关信号匹配技术；实现对甚高频辐射信号中脉冲信息的充分匹配和提取；利用到达时间差技术进行三维定位。本发明可以实现对闪电放电过程超高分辨率的三维定位，定位闪电放电事件的时间分辨率可以达到纳秒量级，空间误差的理论分辨率提升到分米量级。具体步骤如下。
(1)超短基线甚高频天线的多站同步观测布局
本发明将闪电宽带干涉仪(INTF)常用的甚高频(VHF)辐射探测天线共七个，作对称的六边形布局，即以天线A作为中心测站和坐标原点，另有6个相同的天线，依次记为天线B-G，围绕天线A布设，天线B-G中两两相邻的天线的距离为L(称为基线)，天线B-G与天线A的距离也为L。这样，7个天线构成一个高度对称的六边形，参见图1所示。
(2)信号特征分析
将DBM_EEMD(基于集合经验模态分解算法(EEMD)的DBM_EEMD算法，是将待分解信号进行双向双面镜像(double-sided bidirectional mirror，DBM)延拓)算法引入到对闪电产生的甚高频辐射信号的分析中来。分析背景噪声信号、含噪的闪电VHF辐射信号的主要成分。
(3)信号的增采样
对原始信号进行高倍率增采样的预处理，即使用多相滤波器实现对原始信号进行10倍以上(可根据定位研究的目的不同进行调整)的重新采样。对信号的增采样有助于提升信号的时间精度、改变原始信号的不可导特性，提升对信号用DBM_EEMD分解重构的准确性、提升对不同天线上信号做互相关的准确性。
(4)信号的DBM_EEMD降噪
采用DBM_EEMD构造一个带通滤波器，保留探测信号中40-80MHz的信号分量。这样能够有效提高波形匹配的准确性，也有助于获取更为准确的脉冲峰值时间，用于求取不同天线间同一个脉冲信号的时间差，进而实现辐射源的精准定位。
(5)信号的互相关匹配
用广义互相关技术将不同天线上的信号进行匹配，为进一步的脉冲信号识别、匹配做准备。各站的待匹配信号均由主窗口(mainwindow)和两个辅助窗口(auxiliarywindow)构成，中心站的待匹配信号的两个辅助窗口被赋值为0；该技术可以非常有效地提升窗口信号匹配的准确性。以一段192个采样点的信号为例：位于信号中段的长度为64个采样点的主窗口，以及位于主窗口两侧的长度为64个采样点的辅助窗口。具体来说就是以天线A(chA)作为中心站，以每次截取64个采样点作为主窗口，并在主窗口两侧分别构造一个为0值的辅助窗口，而在天线B-G(chX,X＝B…G)上主窗口截取的是与chA同时段的真实信号，不同的是它们的辅助窗口也是向两侧扩展长度为64个采样点的真实信号。辅助窗口的长度取决于INTF天线构成的最长基线的长度，本发明取用64作为辅助窗口的长度。算法中主窗口要保持与辅助窗口相似的权重，因此主窗口的长度也被设定为64个采样点。
(6)微尺度窗口下的脉冲提取
在经过第5步初步匹配后的窗口中，以1个宽度为11ns(可根据信号特征进行调整)的微尺度窗口进行遍历，找出微窗口内存在脉冲峰数量大于等于5个的组合，准确提取脉冲峰的时间。
(7)基于非线性最小二乘算法的辐射源三维定位
包含初始解的求解方法，利用初始解带入非线性最小二乘法得到精确解，可以精确地获取闪电辐射信号发生的三维位置和发生的时刻。
本发明的优势如下。
本发明采用与干涉仪相同参数的辐射信号接收天线，探测到的闪电放电事件的时间分辨率非常高，可以达到纳秒量级，结合信号处理、匹配技术，能够获取闪电放电过程纳秒级分辨率的超精细三维通道和发展特征。
本发明提出的超短基线的天线布局，利用到达时间差定位技术能够精确地获取闪电VHF辐射信号发生的三维空间位置和时刻。
本发明对原始观测信号进行高倍率增采样，进而用DBM_EEMD技术对增采样信号进行优化和带通滤波，进一步提升了VHF辐射信号的时间分辨率以及对波形进行互相关匹配的准确性。
本发明摒弃了以往闪电定位技术中通常采用的质心法，能够实现对闪电VHF辐射信号中逐个脉冲事件的精准定位，即本发明采用的脉冲匹配技术对闪电放电事件的识别定位能力大幅提升。
附图说明
图1为超短基线闪电三维定位系统的天线架设布局。
图2为一段长度为1700个采样点的(a)背景噪声及其(b)频谱，(c)是(b)中水平红框标识部分的局部放大，(d)是(b)中垂直红框的局部放大。
图3为一段长度为1700个采样点的(a)信号及其(b)频谱，(c)是(b)中红框标识部分的局部放大。
图4为信号增采样示意图。
图5为图2a中的背景噪声信号经过带通滤波器后被滤除的部分(Fig.5a)和被保留的部分(Fig.5c)以及它们各自的频谱。
图6为图3a中的闪电VHF辐射信号经过带通滤波器后被滤除的部分(Fig.6a)和被保留的部分(Fig.6c)以及它们各自的频谱。
图7为基于主窗口与辅助窗口相结合的广义互相关波形匹配方法：以chA为中心站，提取时长为352ns的波形(64个采样点)作为chA的主窗口，并在主窗口的两侧分别构造一个同样长度但是值为0的辅助窗口；在天线B-G(chX,X＝B…G)上提取同时段信号用以与chA求取广义互相关系数。
图8为对于经过广义互相关匹配在一起的VHF波形，在微尺度窗口(11ns)内进行脉冲峰的匹配。
图9为天线阵列探测闪电事件示意图(Thomas et al.,2004)。
图10为0.355μs时长窗口内脉冲被成功匹配和定位的情况。其中，(a)为信号经过40-80M带通滤波后的脉冲匹配结果，共有21组脉冲被成功匹配和定位；(b)为信号经过20-80M带通滤波器后的脉冲匹配结果，共有14组脉冲被成功匹配和定位。
具体实施方式
(1)超短基线甚高频天线的多站同步观测布局
本发明将闪电宽带干涉仪(INTF)常用的甚高频(VHF)辐射探测天线按照图1所示布局方案进行架设。该观测方案下以天线A作为中心测站和坐标原点，另有6个相同的天线(编号为B-G)围绕天线A布设，B-G中两两相邻的天线的距离为L，天线B-G与天线A的距离也为L。这样，7个天线构成一个高度对称的六边形。为了确保信号匹配的精度及天线架设的可行性，基线长度L被控制在100米左右(如90-110米，示例中取L＝100m)，是一种超短基线观测系统，该基线长度远远小于国际上现有的三维闪电定位系统。由于采用了干涉仪的探测天线，本发明为了表述方便，依然将所述天线阵列称为INTF阵列。该INTF阵列的7个探测天线性能指标相同，采用的是宽带甚高频(16-88MHz)平板接收天线(天线布局如图1所示)。每个接收器的时间序列波形以180MHz/s的采样速率和16位的采样精度被同步记录，初始时间分辨率为5.5ns，进一步可以升级为500MHz/s的采样速率和16位的采样精度，使信号的初始时间分辨率达到2ns。
此外，与INTF共同组网观测的还有快电场变化测量仪，称快天线(FA)，它以100μs衰减常数测量地面垂直电场的变化.该FA天线在3kHz to>20MHz范围内敏感，与INTF阵列天线以相同的采样率和采样精度被采集卡同步记录，这样的观测方案非常有助于闪电放电物理过程中产生的VHF辐射信号与低频电场波形(通常称为‘sferic’)的精准匹配，特别是对于分析微秒量级放电事件的发展细节具有重要意义。该设备用于闪电放电物理过程的识别，不参与闪电定位。
(2)信号特征分析
类似于Fan et al.(2018,2020)在对低频闪电探测系统定位能力提升中所做的，闪电定位系统定位能力的提升始终是围绕着对探测系统特点和对探测到的电场信号特征的准确分析基础之上。并针对闪电电场信号特点，将经验模态分解(EMD)(与Hilbert变换一起称为Hilbert-Huang Transform，HHT)方法引入到对闪电电场信号(低频/甚低频信号)的分析中来(Fan et al.,2018)，进而提出了基于集合经验模态分解算法(EEMD)的DBM_EEMD算法，即对待分解信号进行双向双面镜像(double-sided bidirectional mirror，DBM)延拓，以优化信号特征并提升算法的降噪性能，特别是能够大幅提升对弱脉冲信号提取的准确性(Fan et al.,2020)。
这里，本发明将DBM_EEMD引入到对闪电产生的甚高频辐射信号的分析中来。如图2、3所示，给出2段INTF探测的VHF波形，作为利用基于DBM_EEMD为内核的HHT进行信号特征分析的示例。首先需要了解的是探测系统获取的背景噪声情况，如图2所示为一段长为1700个采样点(9.44μs)的背景噪声(图2a)及其频谱(图2b)。从图2b中所示的背景噪声的谱带分析可以看到，背景噪声来源为以下几个方面：信号的0线漂移(频率为超低频的噪声信号)；采集频段内的白噪声(图2c,可能叠加了其它来源的弱噪声)；在89MHz附近多个频道的功率较强的广播信号(图2d)。这三类噪声信号对于以广义互相关技术为基础求取时间差进行定位来说会有严重影响。
如图3a所示，为一段同样长为1700个采样点(9.44μs)的VHF辐射信号，整体而言辐射信号较弱(振幅小于±1500，相比于探测信号的最大振幅为±215)。从图3b中的频谱分析中可以看到，背景噪声中的两个强噪声源：0线漂移和广播信号噪声稳定存在，在对探测信号的大量分析对比中也发现背景噪声没有发生大的改变。图3c中是采集频段内的VHF辐射信号的谱分布特征，其中包含覆盖全频段的弱的白噪声信号(可能存在幅值很弱的其它来源噪声)。此外，对比图3c可以看到，在40MHz以下(即探测频段内的相对低频段)存在振幅相对较大的相对低频信号和噪声。在脉冲信号匹配和脉冲峰值时间差的求取过程中，探测频段内相对低频波动对脉冲信号信息的准确提取具有很大干扰。
(3)信号的增采样
由于INTF天线工作在VHF频段，所用的采集卡仅有180M的采样率，采样率与截止频率的比值勉强大于2，即使在未来采用了更高采样率的采集卡，这个比值在短时间内也难以达到10以上，从而探测到的辐射信号如图4中黑色散点所示，是由每5.5纳秒一个散点构成的折线。这种折线波形由于处处不可导，无法用DBM_EEMD直接优化和滤波；这种折线波形由于离散，窄窗口内测量值非常有限，利用互相关方法计算不同天线上的窗口波形的相关系数时误差很大，很难达到预期的信号匹配效果；此外，虽然波形的时间分辨率已经达到了5.5纳秒，但是为了获取分米甚至更高空间分辨率的三维定位结果，这个时间分辨率依然不足。
为了解决以上这些问题，本发明提出对原始信号进行高倍率增采样的预处理方法，即使用多相滤波器实现对原始信号进行10倍以上(可根据定位研究的目的不同进行调整)的重新采样。为方便展示效果，图4中所示为对模拟信号进行10倍增采样的示例，图5及以后的算法实施采用的是对真实闪电VHF辐射信号进行50倍增采样后的信号。
(4)信号的DBM_EEMD降噪
通过以上的信号特征分析，在掌握探测信号的主要特征之后，我们采用DBM_EEMD(Fanet al.,2020)构造一个带通滤波器，只保留探测信号中40-80MHz的信号分量，这样能够有效提高波形匹配的准确性，也有助于获取更为准确的脉冲峰值时间，用于求取不同天线间同一个脉冲信号的时间差，进而实现辐射源的精准定位。
如图5所示为图2a中的背景噪声信号经过带通滤波器后被滤除的部分(图5a)和被保留的部分(图5c)以及它们各自的频谱。图5a中包含了背景噪声中的0线漂移、广播信号和40MHz以下的白噪声分量，从与图2a中原始波形的简单的形态对比可知，用DBM_EEMD构造的带通滤波器能够有效去除噪声的绝大部分分量。带通后的40-80MHz的背景噪声绝对成分已经很小(如图5c所示)，背景噪声的极差(最大值与最小值的差)小于250(采集信号的极差为216)。
对于含噪的VHF探测信号，经过带通滤波器之后，虽然会损失一部分真实信号的分量(如图6a所示)，但是这种对部分信号成分的放弃对于闪电放电事件的准确定位是非常有价值的：信号中仅残余极少量的噪声信号成分，能在最大程度上减小噪声信号的影响；带通滤波后的信号成分相对简单，带宽较窄，可以有效提高信号匹配的准确度；有助于显著提升波形中脉冲信息提取的丰富性和准确性。
(5)信号的互相关匹配
完成对原始信号特征分析并用DBM_EEMD构造带通滤波器对原始信号进行质量控制和信号重构后，本文用广义互相关(Stock et al.,2014)技术将不同天线上的信号进行匹配，为进一步的脉冲信号识别、匹配做准备。
与Stock et al.(2014)的窗口匹配方法截然不同的是，本发明提出和采用的技术路线引入了辅助窗口(auxiliary window)的概念。如图7中所示，一段192个采样点的信号被分割为3个部分：位于信号中段的长度为64个采样点的主窗口(mainwindow)，以及位于主窗口两侧的长度为64个采样点的辅助窗口。具体来说就是以天线A(chA)作为中心站，以每次截取64个采样点作为主窗口，并在主窗口两侧分别构造一个为0值的辅助窗口，而在天线B-G(chX,X＝B…G)上主窗口截取的是与chA同时段的真实信号，不同的是它们的辅助窗口也是向两侧扩展长度为64个采样点的真实信号。辅助窗口的长度取决于INTF天线构成的最长基线的长度，以图1所示布局方案为例，基线长度L＝100m，当闪电产生的辐射信号被INTF天线接收时，产生的时间差不会大于光速传播100m所需的时间，即330ns，180M采样率下(时间分辨率5.5ns)对应大约60个采样点的偏差，为了算法的容错性和通用性，本发明取用64作为辅助窗口的长度。算法中主窗口要保持与辅助窗口相似的权重，因此主窗口的长度也被设定为64个采样点。
这种主窗口与辅助窗口的设定具有多方面的优势。首先，信号不会重复参与定位。这是因为，在中心站(chA)的时间轴上，以64个采样点(352ns)遍历信号，在其它天线的信号中(主窗口与2段辅助窗口)用广义互相关方法找出一段与chA的主窗口信号形态完全一致的信号进行后续的脉冲波形匹配和定位计算，chA上的信号不会被重复利用，因此不会出现被重复定位的信息；另一方面，学者们通常期望通过进一步缩小窗口的大小以获取更多的辐射源信息，但是INTF的基线长度通常限定了同一个放电事件在不同天线间的最大时间差，这也就限定了窗口匹配(Stock et al.,2014)算法所能采用的最小的窗口宽度，而本文提出的这种主窗口与辅助窗口相结合的设定实际上突破了以上限制；更重要的是，在chA上只有主窗口而没有辅助窗口的设定可以显著提升广义互相关进行窗口匹配的准确性，进而提升下一步在更小时间尺度上对脉冲信号进行匹配和信息提取的准确度。
(6)微尺度窗口下的脉冲提取
如图7中所示，chX(X＝B…G)与chA的最大相关系数出现在将chX分别向偏移ΔtAX之后即实现了窗口配对(结果如图8所示)，在以往的干涉仪定位技术Stock et al.(2014)和长基线三维闪电定为系统LMA(Thomas and Ronald,2004)的定位方案中，这个时间差ΔtAX被用于获取“辐射源”方向或位置信息。但是，在以窗口为基础的定位技术中，两个天线上窗口内时间序列的相关性主要取决于一个甚至几个强脉冲信号，并且由于窗口内其它信号的影响，天线间用广义互相关求得的时延通常还会偏离窗口内最强脉冲的峰值时刻。
这样，为了将主窗口内不同天线上的脉冲信号准确地进行配对，还需要更多的步骤和约束条件来实现：
(a)首先在中心站chA的主窗口内找出所有峰值(局部极大值)大于阈值(如图8中的水平虚线所示)的脉冲峰的时刻(以图8中所示chA上的最强脉冲为例，其峰值时刻为Tp)；
(b)以Tp为中心，构建一个宽度为11ns的微尺度窗口同时覆盖chA，chX(X＝B…G)，在chX(X＝B…G)上检测是否有脉冲峰出现在窗口范围内，在11ns的窗口内，当同时有至少5个天线的信号上均有脉冲被检测到时则继续下一步；
(c)对步骤b中初步配对的所有脉冲进行相似性判断，即以配对成功的若干个脉冲各自的峰值时刻(TpA，TpB，…，TpG)为中心，分别截取宽度为11ns的波形，计算脉冲波形之间的相关系数，只有当它们两两间的相关系数均大于0.8以上时，才认为以上所有脉冲信号来自于同一个“放电”事件；
(d)配对成功的脉冲的峰值时刻分别为TpA，TpB，…，TpG，以chA和chB为例，同一个脉冲信号到达A、B天线的时间差为ΔtAB＝TpA-TpB，n个天线构成的天线阵列有
Figure BDA0002946232610000091
组时间差τij，其中n为测站个数，m＝2，当天线个数大于等于5时，能够用到达时间差技术获取脉冲辐射源的三维定位结果。
(7)基于非线性最小二乘算法的辐射源三维定位
在闪电放电事件的定位过程中，需要根据闪电VHF辐射信号到达各个测站的时间差来计算闪电放电事件发生的位置(x，y，z)和时间(t)，共计4个未知参量。在用非线性最小二乘法求解时，要求方程组超解(方程数量大于未知数的个数，即测站数量大于未知数的个数)，因此当有5个及以上数量测站成功匹配到同一个闪电放电事件时，可以计算得到这个脉冲峰到达测站的时间差，可以使用到达时间差法(TOA)计算闪电放电事件的发生位置和时刻，实现三维定位。
如图9所示，在t0时刻，当闪电放电事件发生在(x,y,z)处时，闪电辐射信号以光速c传播，那么，位于(xi,yi，zi)的测站i接收到此事件信号的时间为：
Figure BDA0002946232610000092
对于公式(1)这样的非线性方程，通过联立方程的方式直接求解非常困难，需要将其进行形式转换成为线性方程组再进行求解。
首先，将式(1)进行形式转换，使之成为如下形式：
Figure BDA0002946232610000093
令：
Figure BDA0002946232610000094
r2≡x2+y2+z2 (4)
将以上两式带入到(2)中，可使其转换为以下形式：
Figure BDA0002946232610000095
进一步合并同类项后可以得到：
Figure BDA0002946232610000096
同样可以给出闪电放电信号到达第j个测站的时间方程，将i站与j站方程相减后得可以得到：
Figure BDA0002946232610000101
令：
tij≡ti-tj xij≡xi-xj yij≡yi-yj zij≡zi-zj (8)
使式(7)的形式转换为：
Figure BDA0002946232610000102
进一步定义算子：
Figure BDA0002946232610000103
则有：
xxij+yyij+zzij-c2ttij≡qij (11)
这时，方程(11)是一个关于(x，y,z,t)的线性方程组，它定义了多个四维空间(x,y,z,t)的平面，方程组的解是这些超平面的交汇点。而对于这样的线性方程组而言，关于(x，y，z，t)的解析解是：
Figure BDA0002946232610000104
这里需要注意的是，在实际的观测中，由于站网基线较短，高度项z的系数一般较小，相对于水平位置x和y项的系数通常存在较大的量级差异，且时间项t的系数中包括与光速的乘积，通常是一个量级较大的值，从而导致利用求解超平面方程组的方法得到的高度值z通常会有较大的误差。因而求解超平面方程组的方法给出的解仅可以作为一些求解非线性方程组的优化方法的初始解。
将求解(12)式,获取初始解,将其作为非线性最小二乘迭代计算的初始估计值，这样有利于得到更为合理的数值解析解(x，y，z，t)，并有效提高迭代过程的收敛速度。与国外一些同类闪电定位系统类似，本系统也使用Levenberg-Marquardt算法迭代求解方程组的数值解。对于非线性最小二乘法迭代得到的最优解，以χ2的大小来衡量定位结果的好坏，其中χ2定义为：
Figure BDA0002946232610000105
其中N表示参与定位的测站数量，即有N个测站成功匹配了同一个闪电放电事件，
Figure BDA0002946232610000106
表示脉冲峰值时间提取的误差水平，本发明所述系统时间测量的初始误差为5.5ns，经过增采样和带通滤波后，时间误差水平显著降低，为方便计算和评价，此处将其设定为1ns。
Figure BDA0002946232610000111
表示第i个测站实际接收到闪电放电事件的时间，
Figure BDA0002946232610000112
是由最小二乘迭代计算得到的闪电VHF辐射信号从发生位置(xfit,yfit,zfit)传输到第i个测站的时间。
经过以上步骤，针对本发明提出的超短基线多站VHF辐射信号采集系统获取的数据(1)，在对探测信号特征进行细致分析(2)，利用DBM_EEMD方法对原始信号增采样(3)后进行质量控制(4)；通过主窗口与辅助窗口相结合(5)，利用广义互相关方法实现不同天线信号的波形匹配(6)；在微尺度窗口下(11ns)通过对脉冲信号的阈值约束、相似度约束等实现脉冲信号的准确配对和到达时间差提取(7)；最终用非线性最小二乘法求解得到匹配脉冲辐射源的空间二维坐标(8)，可以得到相比于基于窗口的干涉仪二维定位方法(Stock etal.,2014)和长基线闪电定位系统LMA(Thomas and Ronald,2004)更为丰富的闪电放电信息。以图7中所示长度为64个采样点(时长0.355μs)的chA的主窗口波形为例，基于窗口的定位方法在以32个采样点作为步进步长时，仅可以得到2个辐射源定位结果，当步进步长为64个采样点时，仅能得到1个辐射源定位结果。而采用本文提出的基于全脉冲匹配的定位方法，定位结果的丰富程度则大幅提升，特别是本文引入的DBM_EEMD方法构造的带通滤波器对脉冲信号的提取具有重要作用。如图10a所示，当信号经过40-80M带通滤波后称为相对窄带信号，信号中的脉冲性特征被凸显，0.355μs时长窗口内共有21组符合筛选阈值的脉冲被成功匹配和定位；当信号经过的是20-80M带通滤波器时，信号的上限频率是下限频率的4倍，相对低频信号与相对高频信号的相互叠加，使VHF辐射信号特征更为复杂仅有14组脉冲符合筛选阈值并被成功匹配和定位，这是由于不同频段的信号到达不同天线时，由于到达时间的差不同，与之叠加的信号可能来自不同辐射源或者同一辐射源但是叠加时的相位不同，产生不同的信号增幅或者抵消效果，使信号特征变得更加复杂。信号采集频段外的噪声信号(信号幅度远高于白噪声信号)则对波形匹配、脉冲信号峰值时刻的提取具有严重影响，在不对原始信号进行质量控制的情况下，图10中的示例信号仅有2个脉冲被匹配和定位。

Claims (4)

1.一种基于宽带甚高频辐射信号的超短基线闪电三维定位方法，其特征在于,结合干涉仪的超短基线定位优势和多站到达时间差技术的三维定位能力，采用超短基线多站观测布局；通过对原始信号进行增采样；引入集合经验模态分解技术对增采样信号进行分析和优化；采用基于主窗口与辅助窗口相结合的互相关信号匹配技术,实现对甚高频辐射信号中脉冲信息的充分匹配和提取；利用到达时间差技术进行三维定位；具体步骤为：
(1)超短基线甚高频天线的多站同步观测布局
将闪电宽带干涉仪(INTF)常用的甚高频(VHF)辐射探测天线共七个，作对称的六边形布局，即以天线A作为中心测站和坐标原点，另外6个相同的天线，依次记为天线B-G，围绕天线A布设，天线B-G中两两相邻的天线的距离为L，天线B-G与天线A的距离也为L；这样，7个天线构成一个高度对称的六边形；称L为基线；
(2)信号特征分析
将DBM_EEMD算法用于对闪电产生的甚高频辐射信号进行分析,得到背景噪声信号、含噪的闪电VHF辐射信号的主要成分；这里，DBM_EEMD算法是指基于集合经验模态分解算法(EEMD)，并将待分解信号进行双向双面镜像(DBM)延拓；
(3)信号的增采样
对原始信号进行高倍率增采样的预处理，即使用多相滤波器实现对原始信号进行10倍以上的重新采样，以提升信号的时间精度、改变原始信号的不可导特性，提升对信号用DBM_EEMD分解重构的准确性，提升对不同天线上信号做互相关的准确性；
(4)信号的DBM_EEMD降噪
采用DBM_EEMD构造一个带通滤波器，保留探测信号中40-80MHz的信号分量，以有效提高波形匹配的准确性，获取更为准确的脉冲峰值时间，用于求取不同天线间同一个脉冲信号的时间差，进而实现辐射源的精准定位；
(5)信号的互相关匹配
用广义互相关技术将不同天线上的信号进行匹配，为进一步的脉冲信号识别、匹配做准备；各站的待匹配信号均由主窗口和两个辅助窗口构成，中心站的待匹配信号的两个辅助窗口被赋值为0，以提升窗口信号匹配的准确性；
(6)微尺度窗口下的脉冲提取
在经过步骤(5)初步匹配后的窗口中，以1个宽度为11ns的微尺度窗口进行遍历，找出微窗口内存在脉冲峰数量大于等于5个的组合，以准确提取脉冲峰的时间；
为了将主窗口内不同天线上的脉冲信号准确地进行配对，提取脉冲峰的时间的具体流程为：
(a)首先在中心站天线A的主窗口内找出所有峰值即局部极大值大于阈值的脉冲峰的时刻，记该峰值时刻为Tp；
(b)以Tp为中心，构建一个宽度为11ns的微尺度窗口同时覆盖天线A，天线B，…，天线G，在天线B--天线G上检测是否有脉冲峰出现在窗口范围内，在11ns的窗口内，当同时有至少5个天线的信号上均有脉冲被检测到时，则继续下一步；
(c)对步骤b中初步配对的所有脉冲进行相似性判断，即以配对成功的若干个脉冲各自的峰值时刻TpA，TpB，…，TpG为中心，分别截取宽度为11ns的波形，计算脉冲波形之间的相关系数，只有当它们两两间的相关系数均大于0.8以上时，才认为以上所有脉冲信号来自于同一个“放电”事件；
(d)配对成功的脉冲的峰值时刻分别记为TpA，TpB，…，TpG，对于天线A和天线B，同一个脉冲信号到达A、B天线的时间差为ΔtAB＝TpA-TpB，n个天线构成的天线阵列有
Figure FDA0003738569890000021
组时间差τij，其中n为测站个数，m＝2，当天线个数大于等于5时，用到达时间差技术获取脉冲辐射源的三维定位结果；
(7)基于非线性最小二乘算法的辐射源三维定位
包含初始解的求解，利用初始解带入非线性最小二乘法得到精确解，以精确获取闪电辐射信号发生的三维位置和发生的时刻；具体流程为：
在闪电放电事件的定位过程中，需要根据闪电VHF辐射信号到达各个测站的时间差来计算闪电放电事件发生的位置(x,y,z)和时间(t)，共计4个未知参量；在用非线性最小二乘法求解时，要求方程组超解，因此当有5个及以上数量测站成功匹配到同一个闪电放电事件时，可以计算得到这个脉冲峰到达测站的时间差，使用到达时间差法(TOA)计算闪电放电事件的发生位置和时刻，实现三维定位；
设在t0时刻，当闪电放电事件发生在(x,y,z)处时，闪电辐射信号以光速c传播，那么，位于(xi,yi,zi)的测站i接收到此事件信号的时间为：
Figure FDA0003738569890000022
公式(1)为非线性方程，将其进行形式转换成为线性方程组，再进行求解；
首先，将式(1)进行形式转换，使之成为如下形式：
Figure FDA0003738569890000023
令：
Figure FDA0003738569890000024
r2≡x2+y2+z2 (4)
将以上两式带入到(2)中，使其转换为以下形式：
Figure FDA0003738569890000031
进一步合并同类项后得到：
Figure FDA0003738569890000032
同样可以给出闪电放电信号到达第j个测站的时间方程，将i站与j站方程相减后得得到：
Figure FDA0003738569890000033
令：
tij≡ti-tj xij≡xi-xj yij≡yi-yj zij≡zi-zj (8)
使式(7)的形式转换为：
Figure FDA0003738569890000034
进一步定义算子：
Figure FDA0003738569890000035
则有：
xxij+yyij+zzij-c2ttij≡qij (11)
这时，方程(11)是一个关于(x,y,z,t)的线性方程组，它定义了多个四维空间(x,y,z,t)的平面，方程组的解是这些超平面的交汇点；对于这样的线性方程组，关于(x,y,z,t)的解析解是：
Figure FDA0003738569890000036
为避免有较大的误差，将求解超平面方程组的方法给出的解仅可以作为一些求解非线性方程组的优化方法的初始解；即：
将求解(12)式,获取初始解,将其作为非线性最小二乘迭代计算的初始估计值，这样有利于得到更为合理的数值解析解(x,y,z,t)，并有效提高迭代过程的收敛速度；对于非线性最小二乘法迭代得到的最优解，以χ2的大小来衡量定位结果的好坏，其中χ2定义为：
Figure FDA0003738569890000041
其中，N表示参与定位的测站数量，即有N个测站成功匹配了同一个闪电放电事件，
Figure FDA0003738569890000042
表示脉冲峰值时间提取的误差水平，由于时间测量的初始误差为5.5ns-2ns，经过增采样和带通滤波后，时间误差水平显著降低，为方便计算和评价，此处将其设定为1ns；
Figure FDA0003738569890000043
表示第i个测站实际接收到闪电放电事件的时间，
Figure FDA0003738569890000044
是由最小二乘迭代计算得到的闪电VHF辐射信号从发生位置(xfit,yfit,zfit)传输到第i个测站的时间。
2.根据权利要求1所述的超短基线闪电三维定位方法，其特征在于,步骤(1)中，所述基线L长度控制在100米左右；INTF阵列的7个探测天线性能指标相同，采用宽带甚高频平板接收天线；每个接收器的时间序列波形以180MHz/s～500MHz/s的采样速率和16位的采样精度被同步记录，初始时间分辨率为5.5ns～2ns。
3.根据权利要求2所述的超短基线闪电三维定位方法，其特征在于,与INTF阵列天线共同组网观测的还包括快电场变化测量仪，称快天线(FA)，该快天线(FA)以100μs衰减常数测量地面垂直电场的变化，该FA天线在3kHz到>20MHz范围内敏感，与INTF阵列天线以相同的采样率和采样精度被采集卡同步记录，便于将闪电放电物理过程中产生的VHF辐射信号与低频电场波形的精准匹配。
4.根据权利要求3所述的超短基线闪电三维定位方法，其特征在于,步骤(5)中，所述各站的待匹配信号均由主窗口和两个辅助窗口构成，具体地，将一段192个采样点的信号分割为3个部分：位于信号中段的长度为64个采样点的主窗口，以及位于主窗口两侧的长度为64个采样点的辅助窗口；即以天线A作为中心站，以每次截取64个采样点作为主窗口，并在主窗口两侧分别构造一个为0值的辅助窗口，而在天线B-G上主窗口截取的是与天线A同时段的真实信号，不同的是它们的辅助窗口也是向两侧扩展长度为64个采样点的真实信号；辅助窗口的长度取决于INTF天线构成的最长基线的长度。
CN202110194898.0A 2021-02-21 2021-02-21 基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法 Active CN113075461B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202110194898.0A CN113075461B (zh) 2021-02-21 2021-02-21 基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法
PCT/CN2021/082299 WO2022174486A1 (zh) 2021-02-21 2021-03-23 基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110194898.0A CN113075461B (zh) 2021-02-21 2021-02-21 基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法

Publications (2)

Publication Number Publication Date
CN113075461A CN113075461A (zh) 2021-07-06
CN113075461B true CN113075461B (zh) 2022-08-23

Family

ID=76609410

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110194898.0A Active CN113075461B (zh) 2021-02-21 2021-02-21 基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法

Country Status (2)

Country Link
CN (1) CN113075461B (zh)
WO (1) WO2022174486A1 (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113639970B (zh) * 2021-08-19 2023-11-17 云南电网有限责任公司电力科学研究院 一种评估卫星闪电成像仪探测地面定标能力的方法
CN113820657A (zh) * 2021-10-15 2021-12-21 南通大学 一种闪电三维定位方法
CN114184848B (zh) * 2021-12-03 2023-09-26 中国科学院国家空间科学中心 基于Goertzel算法的星载VHF瞬态信号逐点扫描实时处理方法
CN114428225B (zh) * 2022-04-06 2022-06-14 中国电子科技集团公司第二十九研究所 一种基于准匹配滤波的多辐射源到达角度测量方法及装置
CN114791584B (zh) * 2022-05-20 2023-04-21 南京信息工程大学 基于多站双脉冲间隔差异的信号辐射源定位方法及应用
CN115835467B (zh) * 2022-12-26 2024-03-08 武汉大学 一种基于主被动融合的等离子体三维定位系统及方法
CN116106638B (zh) * 2023-04-11 2023-07-18 中国人民解放军96901部队 一种长基线强电磁脉冲抗干扰定位方法及系统
CN116577568B (zh) * 2023-07-13 2023-10-27 中国气象科学研究院 一种基于哈希值特征的闪电定位方法及装置

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7777676B2 (en) * 2008-07-11 2010-08-17 Toa Systems, Inc. Method of lightning location using time difference of arrival technology with auto correlation
CN102095943B (zh) * 2010-11-25 2013-04-24 中国气象科学研究院 雷电预警方法及装置
CN103235284A (zh) * 2013-03-29 2013-08-07 中国气象科学研究院 一种多站闪电甚高频辐射源三维定位方法和系统
CN104569913B (zh) * 2015-01-30 2017-04-12 武汉大学 一种高精度全闪电定位方法
JP6649160B2 (ja) * 2016-03-31 2020-02-19 株式会社フランクリン・ジャパン 落雷位置標定システム、落雷位置標定装置、及び落雷位置標定方法

Also Published As

Publication number Publication date
WO2022174486A1 (zh) 2022-08-25
CN113075461A (zh) 2021-07-06

Similar Documents

Publication Publication Date Title
CN113075461B (zh) 基于宽带甚高频辐射信号探测的超短基线闪电三维定位方法
Shao et al. Broad band radio interferometry for lightning observations
Wang et al. Beijing Lightning Network (BLNET) and the observation on preliminary breakdown processes
CN105466453B (zh) 一种导航设备在线监测系统及方法
EP0288578A1 (en) Method of detecting object and apparatus therefor
CN104316913A (zh) 多通道接收机实时校准装置及校准与误差补偿方法
Saragiotis et al. Automatic P phase picking using maximum kurtosis and/spl kappa/-statistics criteria
CN101609144B (zh) 闪电辐射源三维定位系统
CN105319487A (zh) 变电站局部放电信号检测和定位系统及方法
KR20140120331A (ko) 부분 방전의 분석 및 위치파악 시스템
Fan et al. Application of ensemble empirical mode decomposition in low-frequency lightning electric field signal analysis and lightning location
CN107329127B (zh) 一种用于雷达系统dbf功能检测的相位线性分析方法及系统
CN106556783A (zh) 一种变电站内基于特高频相控阵原理的局部放电测向方法
CN112965037A (zh) 雷达信号波形不确定度测试系统
CN106324563B (zh) 一种多点无源探测全相位信号分选与到达时间差测量系统
Alammari et al. Lightning mapping: Techniques, challenges, and opportunities
CN112799012B (zh) 一种基于脉冲匹配的宽带干涉仪闪电定位方法及系统
CN109298249B (zh) 一种微波拍波辐射场测量装置及方法
CN107315046A (zh) 一种icp‑ms信号处理方法
CN110531298B (zh) 特高频局部放电传感器性能改进检测方法及系统
CN112505499A (zh) 一种电缆附件绝缘异常的区间划分方法
Mäkelä et al. HF radiation emitted by chaotic leader processes
CN102879098B (zh) 窄带高光谱的快速干涉成像方法及装置
CN105572485A (zh) 一种频谱空间复合显示方法
Fan et al. An improved method for analyzing broadband VHF interferometer lightning observations

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