CN110458777B - 一种基于自适应秩校正的高光谱图像去噪方法、系统及介质 - Google Patents
一种基于自适应秩校正的高光谱图像去噪方法、系统及介质 Download PDFInfo
- Publication number
- CN110458777B CN110458777B CN201910717658.7A CN201910717658A CN110458777B CN 110458777 B CN110458777 B CN 110458777B CN 201910717658 A CN201910717658 A CN 201910717658A CN 110458777 B CN110458777 B CN 110458777B
- Authority
- CN
- China
- Prior art keywords
- hyperspectral image
- clean
- image
- hyperspectral
- adaptive
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 87
- 238000012937 correction Methods 0.000 title claims abstract description 80
- 230000003044 adaptive effect Effects 0.000 title claims description 62
- 239000011159 matrix material Substances 0.000 claims abstract description 71
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 23
- 230000003595 spectral effect Effects 0.000 claims description 21
- 230000000903 blocking effect Effects 0.000 claims description 14
- 239000013598 vector Substances 0.000 claims description 13
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 claims description 11
- 230000003190 augmentative effect Effects 0.000 claims description 10
- 238000012545 processing Methods 0.000 claims description 9
- 238000003860 storage Methods 0.000 claims description 6
- 230000015556 catabolic process Effects 0.000 claims description 5
- 238000004590 computer program Methods 0.000 claims description 5
- 238000006731 degradation reaction Methods 0.000 claims description 5
- 230000008030 elimination Effects 0.000 claims description 5
- 238000003379 elimination reaction Methods 0.000 claims description 5
- 230000003416 augmentation Effects 0.000 claims description 3
- 238000012935 Averaging Methods 0.000 claims description 2
- 238000011084 recovery Methods 0.000 abstract description 10
- 230000006870 function Effects 0.000 description 45
- 238000012360 testing method Methods 0.000 description 9
- 238000001530 Raman microscopy Methods 0.000 description 6
- 238000000354 decomposition reaction Methods 0.000 description 5
- 238000001914 filtration Methods 0.000 description 4
- 230000014509 gene expression Effects 0.000 description 4
- 238000005457 optimization Methods 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 235000008331 Pinus X rigitaeda Nutrition 0.000 description 2
- 235000011613 Pinus brutia Nutrition 0.000 description 2
- 241000018646 Pinus brutia Species 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000002156 mixing Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 230000004083 survival effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 241000208340 Araliaceae Species 0.000 description 1
- 101100235081 Caenorhabditis elegans lem-2 gene Proteins 0.000 description 1
- 101100136092 Drosophila melanogaster peng gene Proteins 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000003795 chemical substances by application Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000008602 contraction Effects 0.000 description 1
- 230000000593 degrading effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/90—Dynamic range modification of images or parts thereof
- G06T5/94—Dynamic range modification of images or parts thereof based on local image properties, e.g. for local contrast enhancement
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于自适应秩校正的高光谱图像去噪方法、系统及介质,本发明高光谱图像去噪方法针对各个矩阵化分块退化高光谱图像、矩阵化分块初始估计的干净高光谱图像建立基于自适应秩校正的高光谱图像去噪模型,通过自适应秩校正惩罚函数实现自适应地抵消核范数对大的奇异值的惩罚,基于l2,1范数消除稀疏噪声,采用交替方向乘子算法进行求解分块干净图像并组合获得最终估计的干净高光谱图像。本发明能够解决现有大多数基于低秩矩阵恢复的高光谱去噪方法利用核范数逼近矩阵秩导致过度惩罚大的奇异值和秩信息固定且需提前定义使得方法不灵活等问题,能够去除高光谱图像中的高斯噪声、脉冲噪声、死线和条带等复杂噪声。
Description
技术领域
本发明涉及高光谱图像去噪技术,具体涉及一种基于自适应秩校正的高光谱图像去噪方法、系统及介质。
背景技术
高光谱图像(HSI)包含丰富的光谱信息,在环境监测、城市规划、地质勘探和农林业生产等领域应用广泛。在成像和传输过程中,因受到多种复杂因素影响,获取的高光谱图像包含各种类型的混合噪声,如高斯噪声、脉冲噪声、条带、死线等。这些复杂的混合噪声不仅影响图像质量,且会对后续应用造成很大困扰。因此,提出一种有效的高光谱图像去噪方法对于高光谱图像处理分析和后续应用意义重大。
传统2D图像去噪方法,如非局部均值、K-SVD、块匹配三维滤波(BM3D)等,采用逐波段处理方式对高光谱图像的每个波段逐一进行去噪。然而,逐波段处理方式忽略了光谱相关性,去噪效果不佳。为解决该问题,研究人员综合利用空间信息和光谱信息进行去噪。卢等人提出一种光谱-空间自适应稀疏表示高光谱图像去噪方法。另外,通过将高光谱图像视为3D张量,采用基于张量分解的方法对高光谱图像进行去噪,如多维维纳滤波与Tucker分解结合方法,PARAFAC方法等。通过充分利用空间相似性和光谱相关性,彭等人提出了一种可分解非局部张量字典学习(TDL)方法。总而言之,上述方法主要用于去除高斯噪声及混合高斯泊松噪声。
但是由于实际中,高光谱图像通常受到高斯噪声、脉冲噪声、死线和条带等多种噪声污染。因而,开发一种有效的高光谱图像去噪方法用以消除多种噪声就显得十分重要。由于干净高光谱图像HSI具有低秩特性,故而研究人员通常将此问题转化为低秩矩阵恢复(LRMR)问题。直接的矩阵秩最小化问题是一个非确定性多项式难(NP-hard)问题,因此现有某些方法采用核范数作为凸代理来逼近矩阵秩,或采用秩逼近矩阵秩。非凸Go分解(GoDec)算法被应用来更好地去除混合噪声。此外,通过利用非独立同分布(non-i.i.d.)混合高斯噪声(NMoGs)假设,结合低秩矩阵分解模型,现有的某些方法提出非独立同分布混合高斯模型,该方法具有去除多种噪声的能力。核范数会过度惩罚大的奇异值,从而使获得的解偏离原始解。而基于秩逼近的低秩矩阵恢复方法,秩信息固定且需要提前定义,这在现实问题中不具有灵活性。
综上所述,现有大多数基于低秩矩阵恢复的高光谱去噪方法利用核范数逼近矩阵秩导致过度惩罚大的奇异值和秩r信息固定且需提前定义使得方法不灵活等问题。
发明内容
本发明要解决的技术问题:为解决现有大多数基于低秩矩阵恢复的高光谱去噪方法利用核范数逼近矩阵秩导致过度惩罚大的奇异值和秩r信息固定且需提前定义使得方法不灵活等问题,提供一种基于自适应秩校正(ARC)的高光谱图像去噪方法、系统及介质,本发明能够去除高光谱图像中的高斯噪声、脉冲噪声、死线和条带等复杂噪声。
为了解决上述技术问题,本发明采用的技术方案为:
一种基于自适应秩校正的高光谱图像去噪方法,实施步骤包括:
3)针对各个矩阵化分块退化高光谱图像Y、矩阵化分块初始估计的干净高光谱图像X建立基于自适应秩校正的高光谱图像去噪模型,所述基于自适应秩校正的高光谱图像去噪模型通过自适应秩校正惩罚函数以实现自适应地抵消核范数对大的奇异值的惩罚,且基于l2,1范数消除稀疏噪声;对所述基于自适应秩校正的高光谱图像去噪模型采用交替方向乘子算法进行求解得到分块干净图像X;
4)将所有的分块干净图像X组合获得最终估计的干净高光谱图像。
可选地,步骤1)的详细步骤包括:
可选地,步骤2)的详细步骤包括:对退化高光谱图像,首先按照步长p以像素(u,v)为中心取大小为q×q×B的子立方体,然后将每个子立方体的每个波段转化为1D列向量来构造矩阵从而得到矩阵化分块退化高光谱图像Y;对初始估计的干净高光谱图像首先按照步长p以像素(u,v)为中心取大小为q×q×B的子立方体,然后将每个子立方体的每个波段转化为1D列向量来构造矩阵从而得到矩阵化分块初始估计的干净高光谱图像X;其中,1≤u≤W,1≤v≤H,其中W、H和B分别表示退化高光谱图像的宽度、高度和光谱带数量,q表示子立方体块大小。
可选地,步骤3)中建立的基于自适应秩校正的高光谱图像去噪模型如下式所示;
上式中,||X||*表示矩阵X核范数,X为分块干净图像,<F(X),X>表示F(X)与X的欧式内积,F(X)表示谱算子,X表示矩阵化分块初始估计的干净高光谱图像,||S||2,1表示矩阵S l2,1范数,λ表示规则化参数,表示矩阵弗罗贝尼乌斯范数,Y表示矩阵化分块退化高光谱图像,S表示矩阵化分块稀疏噪声,σ为与高斯噪声N标准差相关的常数。
可选地,步骤4)中将所有的分块干净图像X组合获得最终估计的干净高光谱图像时,还包括针对分块干净图像X的重叠区域取平均值的方式处理重叠区域像素的步骤。
此外,本发明还提供一种基于自适应秩校正的高光谱图像去噪系统,包括计算机设备,该计算机设备被编程或配置以执行所述基于自适应秩校正的高光谱图像去噪方法的步骤。
此外,本发明还提供一种基于自适应秩校正的高光谱图像去噪系统,包括计算机设备,其特征在于,该计算机设备的存储介质上存储有被编程或配置以执行所述基于自适应秩校正的高光谱图像去噪方法的计算机程序。
此外,本发明还提供一种计算机可读存储介质,该计算机可读存储介质上存储有被编程或配置以执行所述基于自适应秩校正的高光谱图像去噪方法的计算机程序。
此外,本发明还提供一种基于自适应秩校正的高光谱图像去噪系统,包括:
第三程序单元,用于针对各个矩阵化分块退化高光谱图像Y、矩阵化分块初始估计的干净高光谱图像X建立基于自适应秩校正的高光谱图像去噪模型,所述基于自适应秩校正的高光谱图像去噪模型通过自适应秩校正惩罚函数以实现自适应地抵消核范数对大的奇异值的惩罚,且基于l2,1范数消除稀疏噪声;对所述基于自适应秩校正的高光谱图像去噪模型采用交替方向乘子算法进行求解,得到分块干净图像X;
第四程序单元,用于将所有的分块干净图像X组合获得最终估计的干净高光谱图像。
和现有技术相比,本发明具有下述优点:本发明基于自适应秩校正的高光谱图像去噪方法提出的基于自适应秩校正的高光谱图像去噪模型通过自适应秩校正惩罚函数以实现自适应地抵消核范数对大的奇异值的惩罚,自适应秩校正惩罚函数不需要提前定义秩r信息,能够自适应地抵消核范数对大的奇异值的惩罚,从而能对秩进行自适应校正,准确捕捉干净HSI本质低秩信息;本发明基于自适应秩校正的高光谱图像去噪模型通过基于l2,1范数消除稀疏噪声,如死线和条带噪声;因此,本发明能够有效地去除混合噪声,还能有效保留高光谱图像的结构、细节信息,具有非常重要的应用价值。
附图说明
图1为本发明实施例的自适应秩校正的高光谱图像去噪方法流程图。
图2为干净高光谱图像的低秩特性示意图。
图3为本发明实施例中当τ和ε不同取值时的标量函数φ(t)的函数图像。
图4为流行的5种现有去噪方法与本实施例方法的第48个波段去噪结果对比。
图5为流行的5种现有去噪方法与本实施例方法的第57个波段去噪结果对比。
图6为5种现有去噪方法与本实施例方法对频带数为108的去噪结果对比图。
图7为5种现有去噪方法与本实施例方法对频带数为150的去噪结果对比图。
图8为5种现有去噪方法与本实施例方法对频带数为103的去噪结果对比图。
图9为5种现有去噪方法与本实施例方法对频带数为138的去噪结果对比图。
具体实施方式
如图1所示,本实施例基于自适应秩校正的高光谱图像去噪方法的实施步骤包括:
3)针对各个矩阵化分块退化高光谱图像Y、矩阵化分块初始估计的干净高光谱图像X建立基于自适应秩校正的高光谱图像去噪模型,基于自适应秩校正的高光谱图像去噪模型通过自适应秩校正惩罚函数以实现自适应地抵消核范数对大的奇异值的惩罚,且基于l2,1范数消除稀疏噪声;对基于自适应秩校正的高光谱图像去噪模型采用交替方向乘子算法(ADMM)进行求解得到分块干净图像X;
4)将所有的分块干净图像X组合获得最终估计的干净高光谱图像。
步骤1)采用非精确增广拉格朗日乘子算法(IALM)进行初始估计的目的在于得到初始估计的干净高光谱图像;本实施例中,步骤1)的详细步骤包括:
步骤2)的目的在于获取分块高光谱图像的矩阵表示形式,将去噪问题转化为数值优化问题。对于任一退化高光谱图像该退化图像一般可认为是由干净图像、稀疏噪声和高斯噪声三部分构成,即这里与有相同的大小,分别代表干净HSI、稀疏噪声、高斯噪声;W、H和B分别表示HSI的宽度、高度和光谱带数量。在退化高光谱图像中,对于以像素(u,v),(1≤u≤W,1≤v≤H)为中心,大小为q×q×B的子立方体,可通过将子立方体的每个波段转化为1D列向量来构造矩阵在初始估计的干净高光谱图像中,可通过同样的方式在相同位置构造矩阵用于后续秩校正项的构建。相应的高光谱图像退化模型可转换为矩阵形式如式(1)所示:
Y=X+S+N (1)
式(1)中,Y表示矩阵化分块退化高光谱图像,X代表分块干净图像,S和N分别代表矩阵化分块稀疏噪声矩阵和矩阵化分块高斯噪声矩阵,矩阵X,S,N和Y具有相同的大小。
本实施例中,步骤2)的详细步骤包括:对退化高光谱图像首先按照步长p以像素(u,v)为中心取大小为q×q×B的子立方体,然后将每个子立方体的每个波段转化为1D列向量来构造矩阵得到矩阵化分块退化高光谱图像Y;对初始估计的干净高光谱图像首先按照步长p以像素(u,v)为中心取大小为q×q×B的子立方体,然后将每个子立方体的每个波段转化为1D列向量来构造矩阵从而得到矩阵化分块初始估计的干净高光谱图像X;其中,1≤u≤W,1≤v≤H,其中W、H和B分别表示退化高光谱图像的宽度、高度和光谱带数量,q表示子立方体块大小。
步骤3)基于自适应秩校正的高光谱图像去噪模型通过自适应秩校正惩罚函数以实现自适应地抵消核范数对大的奇异值的惩罚,且基于l2,1范数消除稀疏噪声,包括脉冲噪声、死线和条带等的本质稀疏结构。本实施例中,步骤3)中建立的基于自适应秩校正的高光谱图像去噪模型如式(2)所示;
式(2)中,||X||*表示矩阵核范数,X为分块干净图像,<F(X),X>表示F(X)与X的欧式内积,F(·)表示谱算子函数,X表示矩阵化分块初始估计的干净高光谱图像,||S||2,1表示矩阵l2,1范数,λ表示规则化参数,表示矩阵弗罗贝尼乌斯范数,Y表示矩阵化分块退化高光谱图像,S表示矩阵化分块稀疏噪声,σ为与高斯噪声N标准差相关的常数。传统的基于秩-r逼近(rank(X)≤r)的低秩矩阵恢复方法中,秩r信息固定且需要提前定义r的上确界,使得此类在现实问题中不具有灵活性。本实施例基于自适应秩校正的高光谱图像去噪模型利用自适应秩校正惩罚函数||X||*-<F(X),X>来有效逼近矩阵秩,其中,||X||*代表矩阵X的核范数,定义为矩阵奇异值之和,X是对初始估计的干净高光谱图像进行分块操作并矩阵化得到,F(·)为谱算子函数。该惩罚函数能进行自适应秩校正,且能够自适应地抵消核范数对大的奇异值的惩罚。l0范数最小化为NP难问题,将l0范数松弛为l2,1范数来描述稀疏噪声,包括脉冲噪声、死线和条带等的本质稀疏结构,其中l2,1范数为矩阵每一行的l2范数(向量各元素的平方和然后求平方根)之和。l2,1范数对于处理特定结构的稀疏噪声,如死线和条带噪声,非常有效。通过转化,得到基于自适应秩校正的高光谱图像去噪模型式(2)并获得分块后干净图像X,其中λ为一个正则化参数,参见式(2)可知,本实施例基于自适应秩校正的高光谱图像去噪模型的自适应秩校正惩罚函数为核范数减去基于初始估计的线性秩校正项,来有效逼近矩阵秩,该惩罚函数不需要提前定义秩r信息,且能够自适应地抵消核范数对较大奇异值的惩罚,具有自适应秩校正能力。同时,利用l2,1范数来逼近l0范数描述稀疏噪声(包括脉冲噪声、死线和条带等)的本质稀疏结构。利用交替方向乘子算法(ADMM)对得到的凸优化模型进行有效求解。与其它高光谱图像去噪方法相比可知,本实施例所提出的自适应秩校正方法不仅能同时去除高斯噪声、脉冲噪声、死线和条带等复杂噪声,还能有效保留高光谱图像的空间和光谱信息。
图2为干净高光谱图像的低秩先验特性的说明。图2显示干净图像的大部分奇异值都趋近于0,体现了干净图像具有低秩的先验特性。步骤3)中采用提出的自适应秩校正(ARC)方法,对分块高光谱图像去噪,获得分块后干净图像。l0范数最小化问题是NP难问题,且基于秩-r逼近(rank(X)≤r)的低秩矩阵恢复方法,秩r信息固定且需要提前定义r的上确界,使得此类在现实问题中不具有灵活性。利用自适应秩校正惩罚函数来有效逼近矩阵秩,利用l0范数松弛为l2,1范数来描述稀疏噪声,最终得到分块干净图像X。对于任何矩阵约束rank(X)≤r等价于式(3)所示函数表达式:
0=σr+1(X)+…+σn(X)=||X||*-||X||(r) (3)
式(3)中,σr+1(X)~σn(X)分别表示分块干净图像X的第r+1个奇异值~第n个奇异值,||X||*表示矩阵核范数,||X||(r)=σ1(X)+…+σr(X),n=min(q2,B)。
基于当前的迭代值Xt,目标方程(3)获得下一步的迭代值Xt+1,通过解如下凸优化方程如式(4)所示;
式(4)中,Xt+1表示下一步的迭代值,||X||*表示矩阵核范数,Xt表示当前的迭代值,Gt是凸函数||X||(r)在Xt这一点的次梯度。
在实际应用中,预先确定r值是一个挑战性的工作,尤其对那些基于分块的HSI去噪方法来说。另外,对于带噪声的低秩矩阵逼近问题,没有必要通过精确约束秩来得到目标解。因此,本发明提出一种惩罚函数来有效逼近秩函数,自适应进行秩校正。如果初始值与真实矩阵偏差不大,就会在一定程度上包含真实矩阵奇异值或者秩的一些信息。基于这个发现,给定一个与Y具有相同大小的矩阵化分块初始估计的干净高光谱图像X,参见式(2),本实施例提出的自适应秩校正惩罚函数定义为式(5):
||X||*-<F(X),X〉 (5)
式(5)中各符号定义参见式(2)。
F(X)=UDiag(f(σ(X)))VT (6)
式(6)中,F(·)表示谱算子函数,σ(X)表示矩阵X的奇异值,f(·)表示对称函数,Diag用于构造一个对角矩阵,对角线上的元素等于f(σ(X)),U和V的列分别表示奇异值中的左、右奇异向量。
式(7)中,fi(x)表示对称函数f在xi处的函数值,x表示任一向量x,xi表示x的第i个数值,||x||∞表示x的最大范数,φ表示标量函数。其中||x||∞=maxi|xi|,由σ(X)降序排列可得||σ(X)||∞=σ1(X),σ(X)表示矩阵X的奇异值。
式(8)中,φ(t)为标量函数φ,ετ表示ε的τ次方,|t|τ表示t绝对值的τ次方。
如果令F=0,式(5)即为核范数。核范数对每个奇异值进行同等放缩。对图像而言,大的奇异值往往对应更多的结构信息,小的奇异值则对应噪声。因而,我们希望对大的奇异值放缩少或者对小的奇异值值放缩多。给定则有式(9):
式(9)中,左侧各符号的定义参见式(2),右侧各符号中,σi(X)表示X的第i个奇异值,φ(ti)表示函数φ在ti处的函数值。此处n=min(q2,B),q表示矩阵X的行数,B表示矩阵X的列数; 表示矩阵X的第i个奇异值,表示矩阵X奇异值向量的最大范数。
在不同的τ>0和ε>0(τ表示输入参数,ε表示输入参数)的情况下:当自变量t在[0,1]上取值,图3显示了τ、ε在不同取值情况下标量函数φ的各种不同形状。当τ和ε取大于零的不同正数时,针对t∈[0,1],φ(t)值的变化;(a)当ε=0.1,τ取大于零的不同正数时,φ(t)值的变化;(b)当τ=2,ε取大于零的不同正数时,φ(t)值的变化。可以看出,参数τ>0主要控制标量函数φ的形状。对于任何满0<ε<1,如果0<t<ε,标量函数φ的函数值φ(t)趋于0,而如果ε<t≤1,当τ→∞,则标量函数φ的函数值φ(t)趋于1。因此,如果ti接近0,有σi(X)-φ(ti)σi(X)→σi(X)。如果ti接近1,则σi(X)-φ(ti)σi(X)→0。因此,本发明提出的惩罚函数能够自适应地抵消核范数对大的奇异值的惩罚。
此外,本实施例基于自适应秩校正的高光谱图像去噪模型还将l0范数松弛为l2,1范数,来描述稀疏噪声,包括脉冲噪声、死线和条带等的本质稀疏结构。
式(10)中,||S||2,1表示矩阵S的l2,1范数,B表示退化高光谱图像的光谱带数量,q表示子立方体块大小,Sij表示矩阵S第(i,j)个数值。l2,1范数对于处理特定结构的稀疏噪声,如死线和条带噪声,非常有效。基于上述推导过程,可得到步骤3)中建立的基于自适应秩校正的高光谱图像去噪模型如式(2)所示。
步骤3)中对基于自适应秩校正的高光谱图像去噪模型采用交替方向乘子算法(ADMM)进行求解,交替方向乘子算法(ADMM)为现有求解方法,目的是求解目标最小化方程,得到全局最优解。交替方向乘子算法(ADMM)利用目标函数的可分离性,将原问题分解为若干个子问题,然后交替地进行求解。交替方向乘子算法(ADMM)在处理大规模问题和多目标优化问题方面非常有效。基于增广拉格朗日函数的目标函数可表示为:
式(11)中,Λ表示拉格朗日乘子,μ表示惩罚参数,其余各符号定义参见式(2)。目标函数(11)可以通过优化其中一个变量,固定其余变量,迭代求解。特别地,在第t+1迭代,目标函数(11)中的变量更新如下:
Λk+1=Λk+μk(Y-Xk+1-Sk+1) (14)
μk+1=min(t*μk,μmax) (15)
式(12)~(15)中,Xk+1表示分块干净图像X在第k+1次迭代中的值,表示自变量为X,Sk,Λk,μk的增广拉格朗日函数,Sk+1表示矩阵化分块稀疏噪声S在第k+1次迭代中的值,表示自变量为Xk+1,S,Λk,μk的增广拉格朗日函数,Λk+1表示拉格朗日乘子Λ在第k+1次迭代中的值,Λk表示拉格朗日乘子Λ在第k次迭代中的值,μk表示惩罚参数u在第k次迭代中的值,Y表示矩阵化分块退化高光谱图像,μk+1表示惩罚参数u在第k+1次迭代中的值,μmax表示惩罚参数u的最大值,t表示步长,此处步长t取值为1.5。通过简单的推导,对于Xk+1与Sk+1的迭代可表示为式(16)~(17):
式(16)~(17)中,||X||*表示分块干净图像X核范数,F(·)表示谱算子函数,||S||2,1表示矩阵化分块稀疏噪声S的l2,1范数,其余各符号与式(12)~(15)中相同。分别将子问题式(16)~(17)称为X子问题与S子问题。这两个子问题都有封闭解。关于X子问题和S子问题的解可分别由如下【引理1】和【引理2】得出。
【引理1】:给定一矩阵Ek计算如式(18),它的奇异值分解(SVD)可表示为式(19),让αk=1/μk,其奇异值收缩算子服从式(20),且有式(21)和式(22);
式(18)中,Ek表示式(18)右边变量在第k次迭代中的值,Y表示矩阵化分块退化高光谱图像,Sk表示矩阵化分块稀疏噪声S在第k次迭代中的值,F(·)表示谱算子函数,Λk表示朗格朗日乘子Λ在第k次迭代中的值,μk表示惩罚参数u在第k次迭代中的值;
Ek=UkΣk(Vk)T,Σk=Diag(σ(Ek)) (19)
式(19)中,矩阵Ek表示(18)式右边变量在第k次迭代中的值,矩阵Σk对角线上的元素等于Ek的奇异值,Uk和Vk的列分别表示奇异值中的左、右奇异向量,σ(Ek)表示Ek的奇异值,Diag(σ(Ek))用于构造一个对角矩阵,对角线上的元素等于Ek的奇异值。
式(20)中,Xk+1表示分块干净图像X在第k+1次迭代中的值,αk=1/μk,μk表示惩罚参数u在第k次迭代中的值,||X||*表示X的核范数,X表示分块干净图像,Ek表示式(18)右边变量在第k次迭代中的值,||X-Ek||F表示矩阵X-Ek的弗罗贝尼乌斯范数(F-范数)。
式(22)中,Ek表示式(18)右边变量在第k次迭代中的值,αk=1/μk,μk表示惩罚参数u在第k次迭代中的值。
【引理2】:给定一矩阵Wk如式(23),S子问题(17)有最优解,其第j列更新如式(24)。
式(23)中,Y表示矩阵化分块退化高光谱图像,Xk+1表示分块干净图像X在第k+1次迭代中的值,Λk表示拉格朗日乘子Λ在第k+1次迭代中的值,μk表示惩罚参数μ在第k次迭代中的值。
式(24)中,Sk+1(:,j)表示矩阵化分块稀疏噪声S的第j列在第k+1次迭代中的值,||Wk(:,j)||2表示Wk的第j列的2-范数,λ表示规则化参数,μk表示惩罚参数μ在第k次迭代中的值。
步骤3)中对基于自适应秩校正的高光谱图像去噪模型采用交替方向乘子算法(ADMM)进行求解的详细步骤包括:
3.2)使用采用非精确增广拉格朗日乘子算法(IALM)估计获得初始估计图像设置参数值S,Λ,μ,δ=10-4,k=0,其中S表示矩阵化分块稀疏噪声,Λ表示拉格朗日乘子,μ表示惩罚参数,δ表示误差范围,k表示迭代次数;
3.3)判断是否满足收敛,当不收敛时,跳转执行下一步;否则,结束并退出;
3.4)根据式(18)计算矩阵Ek;
3.5)根据式(21)更新Xk+1;
3.6)根据式(23)计算矩阵Wk;
3.7)根据式(24)更新Sk+1;
3.8)根据式(14)更新Λk+1;
3.9)根据式(15)更新μk+1;
3.10)跳转执行步骤3.3)。
步骤3.3)判断是否满足收敛的函数表达式为:
||Y-Xk+1-Sk+1||F≤δ*||Y||F (25)
式(25)中,Y表示矩阵化分块退化高光谱图像,Xk+1表示分块干净图像X在第k+1次迭代中的值,Sk+1表示矩阵化分块稀疏噪声S在第k+1次迭代中的值,||Y-Xk+1-Sk+1||F表示Y-Xk+1-Sk+1的弗罗贝尼乌斯范数(F-范数),δ表示误差范围,||Y||F表示Y的弗罗贝尼乌斯范数(F-范数)。
本实施例中,步骤4)中将所有的分块干净图像X组合获得最终估计的干净高光谱图像时,还包括针对分块干净图像X的重叠区域取平均值的方式处理重叠区域像素的步骤。
下文将通过实验对本实施例基于自适应秩校正的高光谱图像去噪方法进行验证。本实施例中利用Washington DC Mall数据集中的一张波段数为191的256×256的图像进行模拟实验,添加6种不同的噪声情景,详情如下:1)独立同分布高斯噪声:测试数据集的每个频带添加零均值,标准差为0.05的高斯噪声。2)非独立同分布高斯噪声:测试数据集的每个频带添加零均值的高斯噪声。噪声的标准差是0到0.05之间的随机值。3)高斯+脉冲噪声:测试数据集的每个频带添加的高斯噪声如情形2)。另外,随机选取测试集中的30个波段添加脉冲噪声。添加的脉冲噪声比率在0.4到0.5之间随机取值。4)高斯+条带噪声:测试数据集的每个频带添加的高斯噪声如情形2)。随机选取测试集中的30个波段添加条带噪声。条带的数量在20到40之间随机取值。5)高斯+死线噪声:测试数据集的每个频带添加的高斯噪声如情形2)。随机选取测试集中的30个波段添加死线噪声。死线的数量在5到10之间随机取值。6)混合噪声:高斯噪声、脉冲噪声、条带、死线添加到测试数据集的每个频带。这四种噪声的添加方式如情形2)-5)所示。
作为本实施例基于自适应秩校正的高光谱图像去噪方法的对比,采用了5种流行的现有去噪算法,即四维块匹配滤波(Block-Matching 4D filtering,BM4D),低秩矩阵恢复(Low-rank matrix recovery,LRMR),张量字典学习(Tensor dictionary learning,TDL),空间光谱全变差(Spatio-spectral total variation,SSTV),非独立同分布混合高斯(Non-independent and identically distributed mixture of Gaussians,NMoG),与本实施例基于自适应秩校正的高光谱图像去噪方法(ARC)在以上的6种不同的噪声情形下进行比较,具体结果如下表1所示。评价指标有四种,分别是平均峰值信噪比(MPSNR)、平均结构相似度(MSSIM)、平均光谱角距离(MSAD)和计算时间(秒)。
表1:本实施例方法和5种现有去噪方法在四种去噪指标下的效果比较。
从表1可以看出,在独立同分布高斯噪声情景下,本实施例基于自适应秩校正的高光谱图像去噪方法(ARC)略低于TDL;在其他复杂更加符合实际的噪声情形下,本发明所提出的ARC在所有评估指标方面都优于其他方法。
图4为流行的5种现有去噪方法与本实施例方法(ARC)针对Washington DC Mall数据集中在混合噪声情景下的第48个波段去噪结果对比;其中(a)混合噪声情景下的噪声图像,(b)BM4D去噪后的图像,(c)LRMA去噪后的HSI,(d)TDL去噪后的图像,(e)SSTV去噪后的图像,(f)NMoG去噪后的结果,(g)ARC去噪后的图像,及(h)原清洁图像。图5为流行的5种现有去噪方法与本实施例方法(ARC)针对Washington DC Mall数据集中在混合噪声情景下的第57个波段去噪结果对比;其中(a)混合噪声情景下的噪声图像,(b)BM4D去噪后的图像,(c)LRMA去噪后的HSI,(d)TDL去噪后的图像,(e)SSTV去噪后的图像,(f)NMoG去噪后的结果,(g)ARC去噪后的图像,及(h)原清洁图像。图6为5种现有去噪方法与本实施例方法(ARC)在真实的遥感图像数据集AVIRIS Indian Pines中频带数为108的去噪结果对比图;其中(a)被污染的噪声图像,(b)BM4D去噪后的图像,(c)LRMA去噪后的HSI,(d)TDL去噪后的图像,(e)SSTV去噪后的图像,(f)NMoG去噪后的结果,及(g)ARC去噪后的图像。图7为5种现有去噪方法与本实施例方法(ARC)在真实的遥感图像数据集AVIRIS Indian Pines中频带数为150的去噪结果对比图;其中(a)被污染的噪声图像,(b)BM4D去噪后的图像,(c)LRMA去噪后的HSI,(d)TDL去噪后的图像,(e)SSTV去噪后的图像,(f)NMoG去噪后的结果,及(g)ARC去噪后的图像。图8为5种现有去噪方法与本实施例方法(ARC)在真实的遥感图像数据集HYDICE Urban中频带数为103的去噪结果对比图;其中(a)被污染的噪声图像,(b)BM4D去噪后的图像,(c)LRMA去噪后的HSI,(d)TDL去噪后的图像,(e)SSTV去噪后的图像,(f)NMoG去噪后的结果,及(g)ARC去噪后的图像。图9为5种现有去噪方法与本实施例方法(ARC)在真实的遥感图像数据集HYDICE Urban中频带数为138的去噪结果对比图;其中(a)被污染的噪声图像,(b)BM4D去噪后的图像,(c)LRMA去噪后的HSI,(d)TDL去噪后的图像,(e)SSTV去噪后的图像,(f)NMoG去噪后的结果,及(g)ARC去噪后的图像。结合图4~图9也可以看出,本实施例基于自适应秩校正的高光谱图像去噪方法(ARC)可以有效地去除混合噪声并且很好地保留真实的图像结构,具有重要的应用价值。
此外,本实施例还提供一种基于自适应秩校正的高光谱图像去噪系统,包括计算机设备,该计算机设备被编程或配置以执行本实施例前述基于自适应秩校正的高光谱图像去噪方法的步骤。此外,本实施例还提供一种基于自适应秩校正的高光谱图像去噪系统,包括计算机设备,该计算机设备的存储介质上存储有被编程或配置以执行本实施例前述基于自适应秩校正的高光谱图像去噪方法的计算机程序。此外,本实施例还提供一种计算机可读存储介质,该计算机可读存储介质上存储有被编程或配置以执行本实施例前述基于自适应秩校正的高光谱图像去噪方法的计算机程序。此外,本实施例还提供一种基于自适应秩校正的高光谱图像去噪系统,包括:第一程序单元,用于对退化高光谱图像矩阵化后采用非精确增广拉格朗日乘子算法进行初始估计,得到初始估计的干净高光谱图像第二程序单元,用于对退化高光谱图像初始估计的干净高光谱图像分别进行分块和矩阵化处理,分别得到矩阵化分块退化高光谱图像Y、矩阵化分块初始估计的干净高光谱图像X;第三程序单元,用于针对各个矩阵化分块退化高光谱图像Y、矩阵化分块初始估计的干净高光谱图像X建立基于自适应秩校正的高光谱图像去噪模型,基于自适应秩校正的高光谱图像去噪模型通过自适应秩校正惩罚函数以实现自适应地抵消核范数对大的奇异值的惩罚,且基于l2,1范数消除稀疏噪声;对基于自适应秩校正的高光谱图像去噪模型采用交替方向乘子算法进行求解,得到分块干净图像X;第四程序单元,用于将所有的分块干净图像X组合获得最终估计的干净高光谱图像。
以上所述仅是本发明的优选实施方式,本发明的保护范围并不仅局限于上述实施例,凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理前提下的若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。
Claims (6)
1.一种基于自适应秩校正的高光谱图像去噪方法,其特征在于实施步骤包括:
3)针对各个矩阵化分块退化高光谱图像Y、矩阵化分块初始估计的干净高光谱图像建立基于自适应秩校正的高光谱图像去噪模型,所述基于自适应秩校正的高光谱图像去噪模型通过自适应秩校正惩罚函数以实现自适应地抵消核范数对大的奇异值的惩罚,且基于l2,1范数消除稀疏噪声;对所述基于自适应秩校正的高光谱图像去噪模型采用交替方向乘子算法进行求解得到分块干净图像X;所述基于自适应秩校正的高光谱图像去噪模型如下式所示;
上式中,||X||*表示X的核范数,X为分块干净图像,表示与X的欧式内积,表示谱算子,表示矩阵化分块初始估计的干净高光谱图像,||S||2,1表示矩阵l2,1范数,λ表示规则化参数,表示矩阵弗罗贝尼乌斯范数,Y表示矩阵化分块退化高光谱图像,S表示矩阵化分块稀疏噪声,σ为与高斯噪声N标准差相关的常数;
4)将所有的分块干净图像X组合获得最终估计的干净高光谱图像。
4.根据权利要求1所述的基于自适应秩校正的高光谱图像去噪方法,其特征在于,步骤4)中将所有的分块干净图像X组合获得最终估计的干净高光谱图像时,还包括针对分块干净图像X的重叠区域取平均值的方式处理重叠区域像素的步骤。
5.一种基于自适应秩校正的高光谱图像去噪系统,包括计算机设备,其特征在于,该计算机设备被编程或配置以执行权利要求1~4中任意一项所述基于自适应秩校正的高光谱图像去噪方法的步骤。
6.一种计算机可读存储介质,其特征在于,该计算机可读存储介质上存储有被编程或配置以执行权利要求1~4中任意一项所述基于自适应秩校正的高光谱图像去噪方法的计算机程序。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910717658.7A CN110458777B (zh) | 2019-08-05 | 2019-08-05 | 一种基于自适应秩校正的高光谱图像去噪方法、系统及介质 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910717658.7A CN110458777B (zh) | 2019-08-05 | 2019-08-05 | 一种基于自适应秩校正的高光谱图像去噪方法、系统及介质 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110458777A CN110458777A (zh) | 2019-11-15 |
CN110458777B true CN110458777B (zh) | 2021-10-15 |
Family
ID=68484906
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910717658.7A Active CN110458777B (zh) | 2019-08-05 | 2019-08-05 | 一种基于自适应秩校正的高光谱图像去噪方法、系统及介质 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110458777B (zh) |
Families Citing this family (17)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111028172A (zh) * | 2019-12-10 | 2020-04-17 | 浙江工业大学 | 基于无参数的非凸低秩矩阵逼近的高光谱图像去噪方法 |
CN111062888B (zh) * | 2019-12-16 | 2022-02-15 | 武汉大学 | 一种基于多目标低秩稀疏及空谱全变分的高光谱影像去噪方法 |
CN111080555B (zh) * | 2019-12-22 | 2023-08-11 | 北京理工大学 | 一种基于三维拟递归神经网络的高光谱图像降噪方法 |
CN111292266B (zh) * | 2020-01-22 | 2022-07-05 | 武汉大学 | 基于双低秩矩阵分解的gf-5遥感影像混合噪声去除方法 |
CN111397733B (zh) * | 2020-04-23 | 2021-03-02 | 湖南大学 | 一种单/多帧快照式光谱成像方法、系统及介质 |
CN111583149B (zh) * | 2020-05-08 | 2023-02-17 | 西安电子科技大学 | 一种基于l1-l0范数最小化的气动热辐射图像自动校正方法 |
CN111951181B (zh) * | 2020-07-14 | 2024-03-29 | 浙江工业大学 | 基于非局部相似性及加权截断核范数的高光谱图像去噪方法 |
CN112084847B (zh) * | 2020-07-30 | 2024-05-03 | 浙江工业大学 | 基于多通道截断核范数及全变差正则化的高光谱图像去噪方法 |
CN112508807B (zh) * | 2020-11-26 | 2023-07-21 | 电子科技大学 | 一种基于多方向全变分的图像去噪方法 |
CN112700389B (zh) * | 2021-01-13 | 2022-08-19 | 安徽工业大学 | 一种活性污泥微生物彩色显微图像去噪方法 |
CN112785528B (zh) * | 2021-02-01 | 2022-08-02 | 南京信息工程大学 | 一种基于自适应分块旋转滤波的图像去噪方法 |
CN113160069B (zh) * | 2021-02-26 | 2024-02-02 | 桂林电子科技大学 | 一种基于图信号的高光谱图像去噪方法 |
CN113011321B (zh) * | 2021-03-17 | 2022-05-06 | 中南大学 | 一种基于联合字典的光谱信号去噪方法、系统、终端及可读存储介质 |
CN112801881B (zh) * | 2021-04-13 | 2021-06-22 | 湖南大学 | 一种高分辨率高光谱计算成像方法、系统及介质 |
CN113409261B (zh) * | 2021-06-13 | 2024-05-14 | 西北工业大学 | 一种基于空谱特征联合约束的高光谱异常检测方法 |
CN113421198B (zh) * | 2021-06-17 | 2023-10-20 | 南京邮电大学 | 一种基于子空间的非局部低秩张量分解的高光谱图像去噪方法 |
CN116012263A (zh) * | 2023-03-27 | 2023-04-25 | 四川工程职业技术学院 | 一种图像噪声去除方法、装置、存储介质及电子设备 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105069758A (zh) * | 2015-08-21 | 2015-11-18 | 武汉大学 | 一种基于鲁棒低秩张量的高光谱图像去噪方法 |
CN106709881A (zh) * | 2016-12-14 | 2017-05-24 | 上海增容数据科技有限公司 | 一种基于非凸低秩矩阵分解的高光谱图像去噪方法 |
CN108133465A (zh) * | 2017-12-29 | 2018-06-08 | 南京理工大学 | 基于空谱加权tv的非凸低秩松弛的高光谱图像恢复方法 |
US10310074B1 (en) * | 2014-03-27 | 2019-06-04 | Hrl Laboratories, Llc | System and method for denoising synthetic aperture radar (SAR) images via sparse and low-rank (SLR) decomposition and using SAR images to image a complex scene |
-
2019
- 2019-08-05 CN CN201910717658.7A patent/CN110458777B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US10310074B1 (en) * | 2014-03-27 | 2019-06-04 | Hrl Laboratories, Llc | System and method for denoising synthetic aperture radar (SAR) images via sparse and low-rank (SLR) decomposition and using SAR images to image a complex scene |
CN105069758A (zh) * | 2015-08-21 | 2015-11-18 | 武汉大学 | 一种基于鲁棒低秩张量的高光谱图像去噪方法 |
CN106709881A (zh) * | 2016-12-14 | 2017-05-24 | 上海增容数据科技有限公司 | 一种基于非凸低秩矩阵分解的高光谱图像去噪方法 |
CN108133465A (zh) * | 2017-12-29 | 2018-06-08 | 南京理工大学 | 基于空谱加权tv的非凸低秩松弛的高光谱图像恢复方法 |
Non-Patent Citations (2)
Title |
---|
HYPERSPECTRAL IMAGE DENOISING WITH MULTISCALE LOW-RANK MATRIX RECOVERY;Zhihong Huang等;《IGARSS 2017》;20171231;第5442-5445页 * |
基于L21范数正则化矩阵分解的;张怡婷等;《南昌大学学报(理科版)》;20151031;第39卷(第5期);第426-431页 * |
Also Published As
Publication number | Publication date |
---|---|
CN110458777A (zh) | 2019-11-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110458777B (zh) | 一种基于自适应秩校正的高光谱图像去噪方法、系统及介质 | |
CN108133465B (zh) | 基于空谱加权tv的非凸低秩松弛的高光谱图像恢复方法 | |
CN111369487B (zh) | 一种高光谱和多光谱图像融合方法、系统及介质 | |
Guo et al. | Covariance intersection based image fusion technique with application to pansharpening in remote sensing | |
CN105957026B (zh) | 基于非局部相似图像块内部和块间隐性低秩结构的去噪方法 | |
CN109671029B (zh) | 基于伽马范数最小化的图像去噪方法 | |
US20220301114A1 (en) | Noise Reconstruction For Image Denoising | |
CN107680116B (zh) | 一种监测视频图像中运动目标的方法 | |
CN107730482B (zh) | 一种基于区域能量和方差的稀疏融合方法 | |
CN112069919A (zh) | 基于非凸低秩矩阵逼近和全变分正则化的高光谱图像去噪方法 | |
CN110135344B (zh) | 基于加权固定秩表示的红外弱小目标检测方法 | |
CN106447632B (zh) | 一种基于稀疏表示的raw图像去噪方法 | |
CN112598069B (zh) | 基于特征提取和权重系数参数更新的高光谱目标跟踪方法 | |
CN109345563A (zh) | 基于低秩稀疏分解的运动目标检测方法 | |
Chen et al. | Hyperspectral image denoising by total variation-regularized bilinear factorization | |
Yu et al. | Quaternion-based sparse representation of color image | |
CN112598711A (zh) | 一种基于联合光谱降维和特征融合的高光谱目标跟踪方法 | |
CN113421198B (zh) | 一种基于子空间的非局部低秩张量分解的高光谱图像去噪方法 | |
CN109727280B (zh) | 一种基于正交基的高光谱图像丰度估计方法 | |
Malézieux et al. | Dictionary and prior learning with unrolled algorithms for unsupervised inverse problems | |
US20240046602A1 (en) | Hyperspectral image distributed restoration method and system based on graph signal processing and superpixel segmentation | |
CN115731135A (zh) | 基于低秩张量分解和自适应图全变分的高光谱图像去噪方法及系统 | |
CN113160069B (zh) | 一种基于图信号的高光谱图像去噪方法 | |
Bai et al. | Convolutional sparse coding for demosaicking with panchromatic pixels | |
CN111598797A (zh) | 一种基于非凸张量范数和全变分的高光谱图像恢复方法 |
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 |