CN109584256B - 一种基于霍夫直线检测的脉冲星色散值估计方法 - Google Patents

一种基于霍夫直线检测的脉冲星色散值估计方法 Download PDF

Info

Publication number
CN109584256B
CN109584256B CN201811439108.5A CN201811439108A CN109584256B CN 109584256 B CN109584256 B CN 109584256B CN 201811439108 A CN201811439108 A CN 201811439108A CN 109584256 B CN109584256 B CN 109584256B
Authority
CN
China
Prior art keywords
straight line
pulsar
dispersion
hough
theta
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
CN201811439108.5A
Other languages
English (en)
Other versions
CN109584256A (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.)
Beijing Normal University
Original Assignee
Beijing Normal 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 Beijing Normal University filed Critical Beijing Normal University
Priority to CN201811439108.5A priority Critical patent/CN109584256B/zh
Publication of CN109584256A publication Critical patent/CN109584256A/zh
Application granted granted Critical
Publication of CN109584256B publication Critical patent/CN109584256B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/13Edge detection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/20Image enhancement or restoration using local operators
    • G06T5/30Erosion or dilatation, e.g. thinning
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/70Denoising; Smoothing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/60Analysis of geometric attributes
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20048Transform domain processing
    • G06T2207/20061Hough transform

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Geometry (AREA)
  • Image Analysis (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种基于霍夫直线检测的脉冲星DM估计算法,该方法能快速准确地通过对脉冲星观测数据绘制得到的时间‑频率波谱图进行坐标变换,并检测变换后图中的直线,计算出该观测脉冲星的DM值。包括如下步骤:1)通过对时间‑频率二维波谱图进行灰度化、小波阈值去噪、高斯平滑、二值化、腐蚀与膨胀等预处理步骤,得到有明显色散曲线的二值图;2)对二值图中的每一数据点纵坐标频率f变换为新的纵坐标F,变换后得到含有明显直线轨迹的二值图;3)提取二值图边缘,使用霍夫直线变换检测边缘图中的直线;4)由检测出的直线斜率及波谱图表示的时间、频率信息最终计算DM值。本发明能够达到处理步骤简单、DM值计算快速准确的效果。

Description

一种基于霍夫直线检测的脉冲星色散值估计方法
技术领域
本发明属于天文数据处理技术领域,特别涉及一种基于霍夫直线检测的脉冲星DM估计算法。
背景技术
脉冲星信号穿过星际介质到达地球上的射电望远镜接收端的过程中,会与介质中的自由电子发生作用,这种星际介质的色散作用使得同一脉冲高频率的能量比低频率的能量先到达接收端,不同频率f1、f2部分信号到达时间差的公式为:
Figure GDA0003347133400000011
其中DM为色散量,这是一个十分重要的参数,可以用于估计脉冲星的距离和研究星际电子密度的分布。
初步的DM值是在脉冲星搜寻过程中确定的,首先制定尝试消色散过程的DM计划,即选定DM的范围,按一定的步长使用不同的DM值分别进行消色散,得到大量脉冲星的候选体,最后通过人工或者机器学习的方法筛选出真正的脉冲星,同时信噪比S/N峰值对应的DM值为最接近真实值的DM值。整个消色散尝试过程相当于对一个较大取值范围里的DM值进行穷举,是一种简单暴力的方法,此过程耗时多,计算量非常大;更精确的DM值需要通过以高时间分辨率准确测量不同频率的脉冲能量到达射电望远镜的时间延迟,然后通过公式(1)进行计算,缺点在于DM的计算精度受到到达时间差测量精度的限制,同时频率值f1、f2的选择从一定程度上也会影响计算结果,所以通常选择低频率f1和f2,信号能量更强,到达时间差更明显。本发明提出一种将DM值计算转变为求直线斜率的算法,通过对脉冲星观测数据绘制而成的时间-频率二维波谱图使用多种数据预处理、坐标变换、霍夫直线检测等技术,快速准确地计算出DM值。
由于脉冲星本身的信号比较微弱,观测过程又受到大量的RFI干扰和其他噪声的干扰,信号在时间-频率波谱图上的DM色散曲线并不明显,我们使用了很多数据预处理方法来减少背景噪声的干扰,包括小波阈值去噪、高斯滤波、膨胀与腐蚀等。
小波阈值去噪广泛应用于一维信号和二维数据去噪,特别已经在脉冲轮廓的消噪上已经表现出不错的效果。其基本的去噪原理是经小波变换后,信号中的噪声对应幅值较小的系数,且分散在所有系数中,而非噪声信号的能量对应着幅值较大的小波系数且集中在低频部分,通过设置合适的阈值可以去除小波系数较小被认为是噪声的部分,最后进行逆变换重构,实现去噪的目的。其中阈值的选择是影响去噪效果的重要因素,最常用阈值方法是软阈值和硬阈值,两者各自有其优点和缺点,选用哪种阈值需要根据实际的应用场景决定。
高斯滤波是一种常用的线性平滑滤波方法,适用于消除高斯噪声,达到平滑图像的目的。形态学中的膨胀与腐蚀常用于去除二值图中的噪点,通过膨胀与腐蚀的组合实现开运算和闭运算,能有效地去除大量噪点。
在求直线斜率时,我们用到了经典算法霍夫直线检测。霍夫直线检测的结果受噪声干扰较小,具有一定的鲁棒性。它的基本思想是二维数据中的每一点(x,y)在参数空间(ρ,θ)都对应为一条曲线,共线的点所对应的曲线会交于一点,通过统计交于同一点的曲线数目可以得到二维图中对应的直线检测结果。传统的霍夫直线检测计算量较大,很多改进的版本用于改善这一缺点。
发明内容
本发明的目的是解决尝试消色散计算量大、耗时严重的问题。
为此,本发明公开了一种基于霍夫直线检测的的脉冲星DM估计算法,该方法能通过时间-频率二维波谱图快速地计算出对应脉冲星的DM值,包括如下步骤:
1)通过对脉冲星观测数据绘制的含有色散曲线的时间-频率波谱图进行小波软阈值去噪,使得背景噪声有所缓解;
2)对上一步去噪后的时间-频率图进行高斯滤波,二值化,膨胀与腐蚀等进一步数据增强的操作,得到含有脉冲星色散曲线的二值图。
3)对二值图的纵坐标频率f进行坐标变换,令新的纵坐标
Figure GDA0003347133400000021
经过坐标变换后,二值图中的曲线将会成为近似直线。
4)通过边缘检测和霍夫直线检测算法,检测二值图中的直线,得到直线斜率信息。
5)根据波谱图的坐标信息和霍夫直线检测得到的直线斜率,计算此脉冲星的DM值。
附图说明
图1为DM估计算法基本思想。
图2为基于霍夫直线检测的脉冲星DM估计算法流程图。
图3为脉冲星观测数据时间-频率波谱图示例。
图4为小波阈值去噪效果图。
图5为膨胀与腐蚀去噪后的二值图。
图6为经过坐标变换之后的二值图。
图7为边缘检测结果图。
图8为霍夫直线检测的结果图。
具体实施方法
下面结合附图对本发明作进一步说明。
图1为本发明中DM计算方法的示意图。
下面通过图1对该基于霍夫直线检测的脉冲星DM计算方法的基本思想进行详尽的描述。
由公式(1)可知,受星际介质的影响,不同频率部分的脉冲到达接收端的时间延迟与
Figure GDA0003347133400000022
成正比,如果令
Figure GDA0003347133400000023
则有
Figure GDA0003347133400000024
F与t成线性关系。也就是说,如果令时间-频率图中的每一点的频率坐标经过
Figure GDA0003347133400000025
这样的变换,新得到的F-t图中会存在直线形式的脉冲信号轨迹,直线的斜率正好是DM的倒数。由此,脉冲星DM值计算的问题就转变为检测坐标变换后F-t图中的直线进而得到其斜率的问题。
图2为基于霍夫直线检测的脉冲星DM值估计算法的流程图。
下面通过图2详细描述基于霍夫直线检测的脉冲星DM估计算法的各个步骤。
步骤1:由脉冲星观测数据生成某段时间内的时间-频率波谱图,其中每一数据点的颜色深浅表示该点的流量密度。由于脉冲星本身的能量较弱,信号接收过程中又有大量噪声干扰,脉冲星的色散曲线常常淹没在噪声当中,所以所得到的波谱图并没有清晰的曲线轮廓。
步骤2:数据预处理。由步骤1得到的波谱图受噪声的干扰严重,如图3所示,脉冲星信号很弱,与背景并不能明显的分隔开,必须要经过严格的预处理过程才能最终得到良好的检测结果。这次我们选择小波去噪来抑制时间-频率图中的背景噪声,增强色散曲线与背景的对比。使用小波阈值去噪对二维图进行处理,对比了硬阈值去噪和软阈值去噪的实验结果,从图4可以看出硬阈值去噪仍然保留大量的尖锐的噪声,软阈值效果更好。接着,使用高斯滤波平滑波谱图并进行二值化。如图5所示,二值图已有明显的曲线轨迹,但是仍然存在很多细节噪声难以去除。二值图再去噪采用形态学上的膨胀与腐蚀实现开操作和闭操作,使得二值图中的噪点明显减少,余下的少数噪点不足以影响之后的检测结果。
步骤3:坐标变换。我们得到经过预处理的二值图,对各个数据点的纵坐标频率f进行变换。这种坐标变换可以表示为:
I1(t,F)=I2(t,f) (2)
其中,I1为坐标变换后图中的数据点的值,t,F分别为数据点的横纵轴坐标。I2为坐标变换前的时间-频率图,t,f分别为时间、频率坐标,F与f满足
Figure GDA0003347133400000031
图6为坐标变换后的二值图,原始二值图中曲线已经变换为新坐标系下的直线。
步骤4:直线检测。首先进行边缘提取,只对边缘点进行霍夫变换可以大大的减少计算量。由于数据预处理已经去除大量背景噪声,获得的二值图边缘信息较为简单,使用简单的边缘检测算子足以达到要求,如Sobel算子,边缘检测结果如图7。接着对边缘上的每一特征点(t,F)进行霍夫变换。t-F坐标系中的任意一条直线可以用极坐标方程表示:
ρ=tsinθ+Fcosθ (3)
故将每一边缘上特征点(ti,Fi)变换到θ-ρ参数空间中为正弦曲线ρ=tisinθ+Ficosθ,t-F坐标系中共线的特征点交于参数空间中的同一点。实际操作时,θ-ρ参数空间被划分为m×n个单元格(给定θ、ρ的范围,并将两个参数范围分别等分为m、n份)。由于我们检测出的直线对应θ值会用于最终DM值的计算,为得到比较精确的计算结果,这里θ的间隔设置相对较小,比如0.1°、0.2°。特征点变换后的每一条曲线都为其经过的单元格进行投票,投票数最高的前N对(θii),i=1,2,3,…,N,为检测出N条直线对应的参数对。
但是,这里存在一个问题是,应该如何选择其中一个θ值用于最终DM值的计算。当然,可以简单地由多个θ的值求出该脉冲星DM值的范围,这里我们使用D-距离点数目和点平均距离值对检测出的N条直线进行再评测,一般来说,较好的直线周围(这里指点与直线的距离小于D)的信号点数目较多,距离直线的平均距离较小。给定阈值D,对于检测出的N条直线l1,l2,…li…,lN,遍历二值图上的每一个信号点(tj,Fj),有I1(tj,Fj)=1,统计D-距离点数目并计算满足点到直线距离小于D的所有点的平均距离
Figure GDA0003347133400000032
Figure GDA0003347133400000033
其中<(tj,Fj),li>表示点(tj,Fj)到直线li的距离,sum(<(tj,Fj),li>)表示所有满足点到直线li的距离小于D的距离之和,num((tj,Fj))表示满足点到直线距离小于D的点的数目。
以图8霍夫直线检测的结果图为例,霍夫变换后投票数前三的参数对(θ,ρ)为分别(43.3,233),(42.3,236),(45,238),评测结果如表1所示:
表1:直线再评测结果
Figure GDA0003347133400000034
从D-距离点数目和点平均距离结果来看,θ=42.3°的这条直线为更准确的直线.
步骤5:DM计算。由霍夫直线检测得到的θ可以计算出t-F坐标系中直线的斜率为tan(90°-θ),但是波谱图本身可能存在拉伸使得该斜率值并不直接等于脉冲星DM值的倒数。已知波谱图表示的观测时间长度T,观测频率范围f1~f2(f1<f2),二值图的长宽比r,可以计算出实际的直线斜率为
Figure GDA0003347133400000041
如图2的时间-频率图示例,该图为脉冲星J1416-6037,DM值289.2cm-3pc,时长295.6ms的波谱图,观测频率范围为1232MHz~1517MHz,另外
Figure GDA0003347133400000042
此次实验过程中我们将二值图的长宽比统一调整为r=1,则有
Figure GDA0003347133400000043
该估计值与ATNF脉冲星目录中记录值仅差0.1493cm-3pc。
我们用此方法多次处理不同脉冲星的观测数据,同样得到十分接近真实值的DM,证明我们的方法是准确有效的。
总之,本发明的实施例公布的是其较佳的实施方式,但并不限于此。本领域的普通技术人员极易根据上述实施例,领会本发明的精神,并做出不同的引申和变化,但只要不脱离本发明的精神,都在本发明的保护范围之内。

Claims (3)

1.一种基于霍夫直线检测的脉冲星色散量估计方法,其特征在于,包括以下步骤:
1)通过对脉冲星观测数据绘制而成的含有色散曲线的时间-频率波谱图进行小波软阈值去噪,得到背景噪声有所缓解的波谱图;
2)对上一步去噪后的波谱图进行高斯平滑、二值化、膨胀与腐蚀去噪这些进一步数据增强的操作,得到含有脉冲星色散曲线的二值图;
3)对二值图的纵坐标频率f进行坐标变换,令新的纵坐标
Figure FDA0003347133390000011
经过坐标轴变换之后,原二值图中的曲线将会变换为近似的直线;
4)通过边缘提取和霍夫直线检测方法,检测边缘图中的直线,得到直线斜率信息;
5)根据原始波谱图表示的时间、频率坐标信息和霍夫直线检测得到的直线斜率,计算此脉冲星的色散量值;具体包括:由霍夫直线检测得到的θ可以计算出t-F坐标系中直线的斜率为tan(90°-θ),但是波谱图本身可能存在拉伸,横纵轴的单位长度不一致,所以该斜率值并不直接等于脉冲星色散量值的倒数;已知波谱图表示的观测时间长度T,观测频率范围f1-f2;其中f1<f2,和二值图的长宽比值r,可以计算出真实情况下,即横纵坐标单位长度相等时对应的直线斜率,由此可计算出该脉冲星的色散量值为真实直线斜率的倒数。
2.根据权利要求1所述的基于霍夫直线检测的脉冲星色散量估计方法,其特征在于,步骤3)中所述的坐标变换为:
将原始波谱图中纵轴表示的频率f变换为新定义的变量F,且F与f满足关系
Figure FDA0003347133390000012
进行这样的坐标变换的原因在于,由于星际介质的影响,不同频率部分的脉冲星信号到达接受端存在不同的时间延迟,这种时间上的延迟与f-2成正比;通过公式推导,我们很容易得出F与t呈线性关系;所以,经过步骤3)的坐标变换,原波谱图中的色散曲线将会变换为新的F-t图中的直线,直线的斜率大于0,且斜率值正好是色散量的倒数。
3.根据权利要求1所述的基于霍夫直线检测的脉冲星色散量估计方法,其特征在于,步骤4)中所述霍夫直线检测包括对直线检测结果进行再评估的过程,具体如下:
对边缘提取后的边缘图中的每一特征点(ti,Fi)进行霍夫变换;由于t-F坐标系中的任意一条直线用极坐标方程表示为:
ρ=tsinθ+Fcosθ (1)
则每一边缘图上的特征点(ti,Fi)变换到θ-ρ参数空间中为正弦曲线ρ=tisinθ+Ficosθ,t-F坐标系中共线的特征点会交于θ-ρ参数空间中的同一点;实际操作时,θ-ρ参数空间被划分为m×n个单元格;给定θ、ρ的范围,并将两个参数范围分别等分为m、n份;由于我们检测出的直线对应θ值会用于最终色散量值的计算,为得到精确的计算结果,这里θ的间隔设置为0.1°或0.2°;特征点经过变换后的每一条曲线都为其经过的单元格进行投票,投票数最高的前N对(θii),i=1,2,3,…,N,为检测出N条直线对应的参数对;
为了得到精确的色散量估计结果,我们采用D-距离点数目和点平均距离值对检测出的N条直线进行再评测,以得到更准确的直线用于色散量值的计算;评测结果好的直线周围的信号点数目多,距离直线的平均距离小;于是给定阈值D,对于检测出的N条直线l1,l2,…li…,lN,遍历二值图上的每一个信号点(tj,Fj),有I1(tj,Fj)=1,统计D-距离点数目并计算满足点到直线距离小于D的所有点的平均距离
Figure FDA0003347133390000021
最后根据D-距离点数目和平均距离两个评测标准,选定N条直线中最好的直线,将其对应的θ用于步骤5)中的脉冲星色散量计算。
CN201811439108.5A 2018-11-28 2018-11-28 一种基于霍夫直线检测的脉冲星色散值估计方法 Expired - Fee Related CN109584256B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811439108.5A CN109584256B (zh) 2018-11-28 2018-11-28 一种基于霍夫直线检测的脉冲星色散值估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811439108.5A CN109584256B (zh) 2018-11-28 2018-11-28 一种基于霍夫直线检测的脉冲星色散值估计方法

Publications (2)

Publication Number Publication Date
CN109584256A CN109584256A (zh) 2019-04-05
CN109584256B true CN109584256B (zh) 2022-04-12

Family

ID=65925047

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811439108.5A Expired - Fee Related CN109584256B (zh) 2018-11-28 2018-11-28 一种基于霍夫直线检测的脉冲星色散值估计方法

Country Status (1)

Country Link
CN (1) CN109584256B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110427878B (zh) * 2019-07-31 2022-10-14 中国科学院新疆天文台 一种快速射电暴信号识别方法与系统
CN111161169B (zh) * 2019-12-20 2023-06-13 五邑大学 基于霍夫变换的绝对相位噪声去除方法、装置和存储介质
CN111292222B (zh) * 2020-01-22 2023-05-12 中国科学院新疆天文台 一种脉冲星消色散装置及方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106289239A (zh) * 2016-08-15 2017-01-04 中国科学院新疆天文台 一种消除脉冲星到达时间数据中宽频时域干扰的方法
CN108768506A (zh) * 2018-04-08 2018-11-06 四川泰富地面北斗科技股份有限公司 一种基于共同门限的多元多频共视比对授时方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7197381B2 (en) * 2003-12-08 2007-03-27 University Of Maryland Navigational system and method utilizing sources of pulsed celestial radiation

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106289239A (zh) * 2016-08-15 2017-01-04 中国科学院新疆天文台 一种消除脉冲星到达时间数据中宽频时域干扰的方法
CN108768506A (zh) * 2018-04-08 2018-11-06 四川泰富地面北斗科技股份有限公司 一种基于共同门限的多元多频共视比对授时方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A Hierarchical Model with Pseudoinverse Learning Algorithm Optimazation for Pulsar Candidate Selection;Shijia Li et al.;《 2018 IEEE Congress on Evolutionary Computation》;20181004;第1-6页 *
基于深度学习的脉冲星识别系统设计与实现;孙利雷;《中国优秀硕士学位论文全文数据库 基础科学辑》;20170315;第2017年卷(第3期);第A007-3页 *

Also Published As

Publication number Publication date
CN109584256A (zh) 2019-04-05

Similar Documents

Publication Publication Date Title
CN107301661A (zh) 基于边缘点特征的高分辨率遥感图像配准方法
CN109584256B (zh) 一种基于霍夫直线检测的脉冲星色散值估计方法
CN108872962B (zh) 基于分数阶傅里叶变换的激光雷达微弱信号提取和分解方法
CN106296655B (zh) 基于自适应权值和高频阈值的sar图像变化检测方法
CN109919870B (zh) 一种基于bm3d的sar图像相干斑抑制方法
CN108550145B (zh) 一种sar图像质量评估方法和装置
CN106199532B (zh) 基于混合傅立叶-小波分析的探地雷达信号降噪方法
CN102073992B (zh) 一种高分辨率sar卫星图像相干斑去噪方法
CN113325277A (zh) 一种局部放电处理方法
CN103985125B (zh) 一种复杂背景sar图像舰船尾迹检测方法
CN101482617A (zh) 基于非下采样轮廓波的合成孔径雷达图像去噪方法
CN113191979B (zh) 一种分区域sar图像非局部均值去噪方法
CN102222322A (zh) 基于多尺度非局部均值的红外图像背景抑制方法
CN104036461B (zh) 一种基于联合滤波的红外复杂背景抑制方法
CN112327259A (zh) 一种sar图像中干扰信号的消除方法和装置
CN104182941A (zh) 高光谱图像条带噪声去噪方法
CN107784639A (zh) 一种无人机遥感图像改进多边滤波去噪方法
CN103377465B (zh) 基于素描图和核选择的sar图像降斑方法
Gupta et al. A noise robust edge detector for color images using hilbert transform
CN103077507A (zh) 基于Beta算法的多尺度SAR图像降噪方法
Li et al. Noise level estimation method with application to EMD-based signal denoising
Zhang et al. A reverberation noise suppression method of sonar image based on shearlet transform
CN109427042B (zh) 一种提取局部海域沉积层的分层结构和空间分布的方法
CN113689449B (zh) 一种中尺度涡特征参数的反演方法及系统
Zhang et al. Bayesian-based speckle suppression for SAR image using contourlet transform

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: 20220412

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