CN111275652A - 一种去除城市遥感图像中雾霾的方法 - Google Patents

一种去除城市遥感图像中雾霾的方法 Download PDF

Info

Publication number
CN111275652A
CN111275652A CN202010120431.7A CN202010120431A CN111275652A CN 111275652 A CN111275652 A CN 111275652A CN 202010120431 A CN202010120431 A CN 202010120431A CN 111275652 A CN111275652 A CN 111275652A
Authority
CN
China
Prior art keywords
image
scale
remote sensing
color
processing
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.)
Granted
Application number
CN202010120431.7A
Other languages
English (en)
Other versions
CN111275652B (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.)
Xijing University
Original Assignee
Xijing 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 Xijing University filed Critical Xijing University
Priority to CN202010120431.7A priority Critical patent/CN111275652B/zh
Publication of CN111275652A publication Critical patent/CN111275652A/zh
Application granted granted Critical
Publication of CN111275652B publication Critical patent/CN111275652B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10024—Color image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation
    • 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
    • Y02A30/00—Adapting or protecting infrastructure or their operation
    • Y02A30/60—Planning or developing urban green infrastructure

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)

Abstract

一种去除城市遥感图像中雾霾的方法，先输入城市遥感图像，获取城市遥感图像的色彩信息；然后把色彩信息中单波段的灰色图像转换成彩色图像；再对输入的城市遥感图像进行傅里叶变换；然后获取相位特征；再构造对数Gabor滤波器并提取图像幅度特征；然后对幅度特征进行颜色补偿；再获得相位幅度特征；然后用多尺度Retinex理论对城市遥感图像进行雾霾去除处理；再进行局部直方图均衡化处理；接着进行小波多尺度分解和融合处理；最后输出结果；本发明不仅能够有效地去除城市遥感图像中的雾霾，更重要的是能有效保持丰富的细节信息，得到质量明显改善的图像，对比度得到显著提高，具有较高的理论价值和应用前景。

Description

一种去除城市遥感图像中雾霾的方法
技术领域
本发明属于遥感技术和图像处理技术及计算机视觉技术领域，具体涉及一种去除城市遥感图像中雾霾的方法。
背景技术
雾霾已成为一种常见的恶劣天气现象，其对户外监测和遥感探测***带来了较大的挑战和困难，对遥感图像进行雾霾的去除处理，减少其影响，提高图像质量，已是遥感图像预处理的重要基础内容。
图像雾霾去除处理起源于室外图像预处理，目前发展了许多行之有效的理论和方法。这些方法均是从不同的层次对图像中的雾霾进行滤除处理，大体上可以分成四类：第一类方法是基于滤波的图像增强方法，如直方图均衡法，Retinex理论和双边滤波；第二类方法是基于大气物理模型的图像恢复算法，典型方法有暗通道先验法；第三类方法是基于融合理论的雾霾去除算法；第四类方法是基于各种机器学习的算法。
直方图是一种经典的图像灰度值范围拉伸方法，但是图像雾霾去除的效果不显著；整体上Retinex理论能够很好地去除图像中的雾霾，但处理的结果在色彩上存在偏灰和失真的现象，并且对于增强雾霾图像中较亮处的细节信息有一定的困难；利用图像增强滤波理论去除雾霾，没有考虑雾霾图像质量降低的原因。
基于大气散射物理模型的雾霾去除是在基于这样的两个假设前提下，即晴天自然清晰图像的对比度高于雾霾笼罩下的图像对比度，以及环境光照量取决于物体在图像中深度的变化，而且还需要估计环境大气光量参数和大气透射率参数，参数估计的精度直接影响雾霾去除效果。
融合策略处理图像，关键技术是融合规则的确定，不同规则产生不同结果。
以深度卷积神经网络技术为代表的深度学习理论在图像分类和目标检测中得到了成功应用，最大困难是图像样本的获取与训练；由于要产生大量的训练样本特征子图像集，所以整个处理过程比较复杂，而且速度缓慢，实时性差。
遥感图像雾霾去除既具有挑战性，又具有非常重要的应用价值，逐渐得到了许多学者的重视。虽然图像雾霾去除算法很多，而且均能获得好的处理效果，但是处理的对象通常都是室外彩色图像或者RGB图像。如果直接用这些方法处理航天或航空遥感图像，并非能取得满意的结果。因为遥感图像与普通室外彩色图像存在较大差异。第一，他们的分辨率差别非常大，遥感图像分辨率低，室外图像分辨率高；第二，室外图像是彩色图像，而遥感图像有灰度、彩色或多光谱图像；第三，有部分室外图像包含有天空背景，而且景深比较长，即景距跨越比较大，因此室外图像的处理重点在远景部分；第四遥感图像包含的是地球表面场景的灰度分布情况，景深比较小。虽然雾霾能够影响图像的质量和清晰度，但是它改变不了遥感图像中地物目标的灰度分布规律，也就是说，不管雾霾存在不存在，是否给遥感图像带来影响，遥感图像中地物目标的相位信息是不变的。
发明内容
为了克服上述现有技术的缺点，本发明的目的在于提供了一种去除城市遥感图像中雾霾的方法，不仅能够有效地去除城市遥感图像中的雾霾，更重要的是能有效保持丰富的细节信息，得到质量明显改善的图像，对比度得到显著提高，具有较高的理论价值和应用前景。
为了达到上述目的，本发明采取的技术方案为：
一种去除城市遥感图像中雾霾的方法，包括以下步骤：
步骤1：输入城市遥感图像；
步骤2：获取城市遥感图像的色彩信息；
步骤3：把色彩信息中单波段的灰色图像转换成彩色图像；
步骤4：对输入的城市遥感图像进行傅里叶变换；
步骤5：获取相位特征；
步骤6：构造对数Gabor滤波器并提取图像幅度特征；
步骤7：对幅度特征进行颜色补偿；
步骤8：获得相位幅度特征；
步骤9：用多尺度Retinex理论对城市遥感图像进行雾霾去除处理；
步骤10：进行局部直方图均衡化处理；
步骤11：进行小波多尺度分解和融合处理；
步骤12：输出结果。
所述的步骤1中输入的遥感图像是单波段灰色图像，或是多光谱图像或多个波段合成的彩色图像，或是由CCD相机获得的RGB真彩色图像，它们都是在雾霾存在的条件下获得的城市区域遥感图像。
所述的步骤2中获取图像的色彩信息的具体步骤为：
步骤2.1：判断输入的图像是彩色图像还是全色的灰度图像，并保存相关的色彩信息，为后续色彩补充和输出结果提供信息；
步骤2.2：输入的是灰度图像，直接进入下一个处理步骤，即步骤3；
输入的是多光谱彩色图像，跳过步骤3，直接进入步骤4。
所述的步骤3中把单波段的灰色图像转换成彩色图像的具体步骤为；
步骤3.1：对输入的全色灰度遥感图像进行彩色增强处理，输入的单波段灰度遥感图像用I(m,n)表示，转换后的RGB假彩色遥感图像用IC(m,n)表示，转换的过程如式(1)所示，就是把灰度图像I(m,n)的各个像素值同时分别赋予给假彩色RGB图像IC(m,n)的三个通道，即ICR(m,n)通道，ICG(m,n)通道和ICB(m,n)通道；
Figure BDA0002392797210000051
所述的步骤4中对输入的城市遥感图像进行傅里叶变换，要对原始遥感图像进行二维傅里叶变换，把图像从空间域转换到频率域。
所述的步骤5中获取相位特征，由于遥感图像经傅里叶变换后，在频率域中变成了相应的复数，包含实部和虚部，以及每个像素的相位，直接提取相位角度信息，但要进行换算，把度转换成弧度；获得的相频特征将与后面获得的幅频特征进行卷积和融合处理，以获得最后的相位幅度特征。
所述的步骤6中构造对数Gabor滤波器并提取图像幅度特征的具体步骤为：
步骤6.1：在频率域直接获取每个像素的幅度值，遥感图像幅度特征的提取与滤波是同时进行，在频率域中，通过构造对数Gabor滤波器对频域中的幅度特征进行滤波，把获得的结果进行傅里叶逆变换，就获得了遥感图像的幅度特征；
步骤6.2：构造Gabor滤波器，对Gabor滤波器设计5个尺度，每个尺度设计6个方向，总共有30个滤波器，把幅度特征经滤波器处理，获得30幅特征图；
步骤6.3：把获得的30幅特征图与步骤5获得的相位信息分别进行卷积运算处理，然后把不同尺度不同方向的处理结果进行加权求和，这里设置各权值相等，即进行求和后再取均值运算，结果作为最终的幅度特征。
所述的步骤7中对幅度特征进行颜色补偿的具体步骤为：
利用步骤2获取的颜色信息对幅度特征进行色彩补偿处理，处理的模型公式如下：
Figure BDA0002392797210000061
Figure BDA0002392797210000062
式中FAi(m,n)和
Figure BDA0002392797210000063
分别表示补偿前后的幅度特征，i表示颜色通道数，Ci(m,n)表示第i个颜色通道的色彩补偿系数，β是一个补偿常数；在其他条件不变的情况下，Ci(m,n)值越大，图像颜色失真越大，同时图像偏暗；当其值较小时，颜色失真小，图像偏亮；Ci(m,n)＝20。
所述的步骤8中获得相位幅度特征的具体步骤为：
按式(4)把相位特征和幅度特征进行融合处理，以获得最终的相位幅度特征图，
FAP(m,n)＝α·FA(m,n)+β·FP(m,n) (4)
式中FAP(m,n)表示融合后的幅度相位特征，FA(m,n)表示幅度特征，FP(m,n)表示相位特征，α和β表示加权系数，分别设为0.45和0.55。
所述的步骤9中用多尺度Retinex理论对城市遥感图像进行雾霾去除处理的具体步骤为：
步骤9.1：获取的图像I(m,n)看作是环境入射光和物体反射相乘的结果，那么它可以用式(5)来表示。
I(m,n)＝L(m,n)·R(m,n) (5)
式中I(m,n)表示用光学传感器获取的图像，即受到了雾霾影响而降质的图像，(m,n)表示像素点的空间位置；L(m,n)表示入射的环境光，即要消除的干扰信息；R(m,n)表示物体表面反射光，即需保留的图像信息；
对公式(5)两边取对数，并进行调整，得到式(6)，
r(m,n)＝log[R(m,n)]＝log[I(m,n)]-log[L(m,n)] (6)
其中L(m,n)用下式进行计算，
L(m,n)＝G(m,n)*I(m,n) (7)
式中的G(m,n)为高斯函数，其定义如下：
Figure BDA0002392797210000071
式中参数σ既为高斯函数的标准偏差，又为其尺度参数；
多尺度Retinex(MSR)算法就是把不同单个尺度Retinex(SSR)算法处理结果进行加权求和，然后把它作为最后的处理结果，其定义如式(9)所示，
Figure BDA0002392797210000072
式中的i为颜色通道数，i∈{1,2,3}，对应于彩色图像的R通道、G通道和B通道；M表示尺度参数值的个数，j表示参数值M中的第j个尺度参数σ，rji(m,n)表示第j尺度下第i颜色通道的反射分量，ri(m,n)表示第i颜色通道经Retinex算法处理后的反射分量，wj为尺度参数值的权值；
步骤9.2：MSR算法中尺度参数σ的个数M及其数值的设计，设置了两种多尺度参数个数和数值方案，具体如下：
MSR算法中尺度参数σ的第一种设计方案：此方案中只设置了6个尺度参数，第一尺度20、第二尺度40、第三尺度60、第四尺度80、第五尺度100、第六尺度120；
MSR算法中尺度参数σ的第二种设计方案：此方案中只设置了9个尺度参数，第一尺度32、第二尺度64、第三尺度128、第四尺度256、第五尺度512、第六尺度1024、第七尺度2048、第八尺度4096、第九尺度8192；
两种方案结果的融合：把两种参数设置方案分别执行处理，然后把各自获得的处理结果进行简单的融合处理，即求和取均值，然后把其作为最终结果。
所述的步骤10中进行局部直方图均衡化处理的具体步骤为：
用MSR算法对图像进行雾霾去除，处理后整体色调偏暗；接着用直方图均衡化进行处理；用直方图原理对图像进行处理时，先把RGB模型图像转换成HSI模型图像，然后对HSI模型中的I分量进行处理；处理结束后，再把HIS模型图像转换成RGB图像；
RGB图像和HSI模型图像转换过程如下：
由RGB图像获取HSI图像的H、S和I分量图的计算公式如下
Figure BDA0002392797210000091
其中θ按式(11)计算，
Figure BDA0002392797210000092
Figure BDA0002392797210000093
Figure BDA0002392797210000094
由HSI图像转换成RBG图像的具体过程如下：
如果0°≤H＜120°，则R、G、B分量由式(14)获取；
Figure BDA0002392797210000095
如果120°≤H＜240°，则R、G、B分量由式(15)进行计算，
Figure BDA0002392797210000096
如果240°≤H＜360°，则R、G、B分量由式(16)进行计算，
Figure BDA0002392797210000097
所述的步骤11中进行小波多尺度分解和融合处理的具体步骤为：
步骤11.1：采用二维平稳小波对图像进行融合处理；
步骤11.2：低频系数的融合，用二维平稳小波分别对相位幅度特征图和滤波结果图进行多尺度分解，在各个分解尺度上再分别进行融合处理，把融合后的结果进行小波逆变换，就获得了最后的结果；采用像素级进行融合处理策略，两幅图像被分解后，分别获得了若干个高频系数和低频系数子图像；对于低频系数子图像，采取系数加权的形式进行融合处理；对于高频系数子图像采取抓大舍小的原则进行处理；设I1(m,n)和I2(m,n)分别表示两幅按不同方法处理后的遥感图像，IF(m,n)表示融合后的图像，低频系数的加权融合规则数学模型如式(17)所示，
IF(m,n)＝w1·I1(m,n)+w2·I2(m,n) (17)
式中的w1和w2分别表示图像I1(m,n)和图像I2(m,n)的融合权值，即加权系数，而且w1+w2＝1；如果w1＝w2＝0.5，则是平均加权，即均值处理；
步骤11.3：高频系数的融合，对于每个分解尺度上高频系数子图像的融合采用取大放小的原则进行融合处理，其融合规则数学模型如式(18)所示，
Figure BDA0002392797210000101
所述的步骤12中输出结果的具体为：根据获取的图像色彩信息，如果最初输入的是灰度图像，还需转换成灰度图像输出；如果是彩色图像，则直接输出相应的结果。
本发明与现有技术相比的优点在于：
本发明是一种基于遥感图像相位一致性、图像统计特性和多尺度理论的城市遥感图像雾霾去除方法，具体的优点如下：
(1)针对性强：本发明针对城市遥感图像提出的发明，能有效去除城市遥感图像中的雾霾，不论是但波段灰度图像还是多波段的彩色图像，都能取得很好的效果，而且能获得增强的丰富的细节和对比度信息。如果图像中包含大量植被信息、天空区域和海洋区域，效果就不理想，那不是本发明的应用范围，因此，本发明的针对性和应用性都很强。
(2)设计思路新：本发明通过获取图像的相位一致性来增强图像的细节和对比度，通过多尺度Retinex理论和直方图均衡化级联处理，达到有效去除雾霾的目的，同时能尽可能保留原始图像的色调不失真，最后经二维平稳小波多尺度融合处理，获得的结果更精确。
(3)处理效果好：对不同的城市遥感图像和不同的处理方法进行了比较，不论从视觉效果还是客观评价参数，表明本发明处理的结果非常不错，既能有效去除雾霾和光晕现象，又能保持丰富的几何细节信息，同时具有高的对比度和清晰度。
附图说明
图1为本发明方法的流程图。
图2为本发明方法与其他方法的实验结果图，其中A-J分别代表不同的城市遥感原始雾霾图像，a-e分别代表不同方法获得的实验结果；a.暗通道先验法(DCP)方法,b.亮度保持模糊直方图(BPDFHE)方法，c.直方图均衡化(HE)方法，d.多尺度模型及直方图特性(MSMHC)方法,e.本发明方法。
具体实施方式
下面结合实施例和附图对本发明做详细描述。
参见图1，一种去除城市遥感图像中雾霾的方法，包括以下步骤：
步骤1：输入城市遥感图像；
步骤2：获取城市遥感图像的色彩信息；
步骤3：把色彩信息中单波段的灰色图像转换成彩色图像；
步骤4：对输入的城市遥感图像进行傅里叶变换；
步骤5：获取相位特征；
步骤6：构造对数Gabor滤波器并提取图像幅度特征；
步骤7：对幅度特征进行颜色补偿；
步骤8：获得相位幅度特征；
步骤9：用多尺度Retinex理论对城市遥感图像进行雾霾去除处理；
步骤10：进行局部直方图均衡化处理；
步骤11：进行小波多尺度分解和融合处理；
步骤12：输出结果。
所述的步骤1中输入的遥感图像是单波段灰色图像，或是多光谱图像或多个波段合成的彩色图像，或是由CCD相机获得的RGB真彩色图像，它们都是在雾霾存在的条件下获得的城市区域遥感图像，所以这些输入的遥感图像质量比较低，对比度低，视觉效果差。
所述的步骤2中获取图像的色彩信息的具体步骤为：
步骤2.1：判断输入的图像是彩色图像还是全色的灰度图像，并保存相关的色彩信息，为后续色彩补充和输出结果提供信息；
步骤2.2：输入的是灰度图像，直接进入下一个处理步骤，即步骤3；
输入的是多光谱彩色图像，跳过步骤3，直接进入步骤4。
所述的步骤3中把单波段的灰色图像转换成彩色图像的具体步骤为；
步骤3.1：为了便于和多波段的彩色遥感图像进行相同的处理步骤和过程，对输入的全色灰度遥感图像进行彩色增强处理，输入的单波段灰度遥感图像用I(m,n)表示，转换后的RGB假彩色遥感图像用IC(m,n)表示，转换的过程如式(1)所示，就是把灰度图像I(m,n)的各个像素值同时分别赋予给假彩色RGB图像IC(m,n)的三个通道，即ICR(m,n)通道，ICG(m,n)通道和ICB(m,n)通道；
Figure BDA0002392797210000131
像暗通道先验法这样的典型雾霾去除算法，要求输入的图像不仅是彩色图像，而且三个通道图像的像素值要求不相等，这是基于大气散射模型的雾霾去除算法的不足。本发明根据实际的遥感图像，通过式(1)把灰度图像转换成彩色图像，就可以直接进行处理了。本发明不仅能处理彩色多光谱图像，同样也能处理单波段的灰度图像，表明本发明具有较宽的适用范围和较好的普适性。
所述的步骤4中对输入的城市遥感图像进行傅里叶变换，充分利用遥感图像中地物目标固有的相位信息，相位信息是相位一致性特征的重要内容，它的获取是在频域中完成的，因此要对原始遥感图像进行二维傅里叶变换，把图像从空间域转换到频率域。
所述的步骤5中获取相位特征，相位一致性特征是城市遥感图像的重要内容和内涵，它包含了非常丰富而且显著的边缘和轮廓信息，对于城市建筑和目标检测具有非常重要的意义，遥感图像的相位一致性特征提取跟步骤6相关。由于遥感图像经傅里叶变换后，在频率域中变成了相应的复数，包含实部和虚部，以及每个像素的相位，直接提取相位角度信息，但要进行换算，把度转换成弧度；获得的相频性特征将与后面获得的幅频特征进行卷积和融合处理，以获得最后的相位幅度特征。
所述的步骤6中构造对数Gabor滤波器并提取图像幅度特征的具体步骤为：
步骤6.1：在频率域直接获取每个像素的幅度值，遥感图像幅度特征的提取与滤波是同时进行，在频率域中，通过构造对数Gabor滤波器对频域中的幅度特征进行滤波，把获得的结果进行傅里叶逆变换，就获得了遥感图像的幅度特征；
步骤6.2：构造Gabor滤波器，对Gabor滤波器设计5个尺度，每个尺度设计6个方向，总共有30个滤波器，把幅度特征经滤波器处理，获得30幅特征图；
步骤6.3：把获得的30幅特征图与步骤5获得的相位信息分别进行卷积运算处理，然后把不同尺度不同方向的处理结果进行加权求和，这里设置各权值相等，即进行求和后再取均值运算，结果作为最终的幅度特征，这样将使特征更加清晰。
所述的步骤7中对幅度特征进行颜色补偿的具体步骤为：
幅度特征经过频域的滤波处理和相位信息的卷积处理后，会产生一定程度的色彩失真，利用步骤2获取的颜色信息对幅度特征进行色彩补偿处理，处理的模型公式如下：
Figure BDA0002392797210000151
Figure BDA0002392797210000152
式中FAi(m,n)和
Figure BDA0002392797210000153
分别表示补偿前后的幅度特征，i表示颜色通道数，Ci(m,n)表示第i个颜色通道的色彩补偿系数，β是一个补偿常数；补偿系数Ci(m,n)是一个重要参数，其值对处理结果的色彩有较大的影响。在其他条件不变的情况下，Ci(m,n)值越大，图像颜色失真越大，同时图像偏暗；当其值较小时，颜色失真小，图像偏亮；对不同的景像，其影响效果程度也不一样。经过大量的实验，并对实验结果进行比较分析，Ci(m,n)∈[10,30]比较合适，对不同的景像图像均能获得比较好的结果；在本发明中，把其值设置为20，即Ci(m,n)＝20。
所述的步骤8中获得相位幅度特征的具体步骤为：
按式(4)把相位特征和幅度特征进行融合处理，以获得最终的相位幅度特征图，
FAP(m,n)＝α·FA(m,n)+β·FP(m,n) (4)
式中FAP(m,n)表示融合后的幅度相位特征，FA(m,n)表示幅度特征，FP(m,n)表示相位特征，α和β表示加权系数，分别设为0.45和0.55。
所述的步骤9中用多尺度Retinex理论对城市遥感图像进行雾霾去除处理的具体步骤为：
步骤9.1：获取的图像I(m,n)看作是环境入射光和物体反射相乘的结果，那么它可以用式(5)来表示。
I(m,n)＝L(m,n)·R(m,n) (5)
式中I(m,n)表示用光学传感器获取的图像，即受到了雾霾影响而降质的图像，(m,n)表示像素点的空间位置；L(m,n)表示入射的环境光，即要消除的干扰信息；R(m,n)表示物体表面反射光，即需保留的图像信息；
对公式(5)两边取对数，并进行调整，得到式(6)，
r(m,n)＝log[R(m,n)]＝log[I(m,n)]-log[L(m,n)] (6)
其中L(m,n)用下式进行计算，
L(m,n)＝G(m,n)*I(m,n) (7)
式中的G(m,n)为高斯函数，其定义如下：
Figure BDA0002392797210000171
式中参数σ既为高斯函数的标准偏差，又为其尺度参数；
式(8)中的尺度参数σ取不同值，对处理效果有较大的影响，而且通过对单个参数值进行调整和选择，很难使细节保留和颜色保真同时达到最佳状态，这就是单尺度Retinex(SSR)算法的不足；为了克服SSR算法的缺陷，达到既能去除雾霾，又能保护细节和色彩，可以通过设置M个参数σ值来完成，即多尺度Retinex(MSR)算法。MSR算法就是把不同单个尺度参数的SSR算法处理结果进行加权求和，然后把它作为最后的处理结果，其定义如式(9)所示，
Figure BDA0002392797210000172
式中的i为颜色通道数，i∈{1,2,3}，对应于彩色图像的R通道、G通道和B通道；M表示尺度参数值的个数，j表示参数值M中的第j个尺度参数σ，rji(m,n)表示第j尺度下第i颜色通道的反射分量，ri(m,n)表示第i颜色通道经Retinex算法处理后的反射分量，wj为尺度参数值的权值；
步骤9.2：MSR算法中关键技术是尺度参数σ的个数M及其数值的设计，不同数目的参数个数和不同取值，都会对滤波效果产生影响。本发明中设置了两种多尺度参数个数和数值方案，具体如下：
MSR算法中尺度参数σ的第一种设计方案：为了保持更多的细节，此方案中只设置了6个尺度参数，具体情况如表1所示，它们的数值比较偏小；这是因为当尺度参数值取较小值时，有利于保护更多细节信息。但是，同时会产生雾霾去除效果不佳及出现光晕和色彩失真的现象；
表1 6个尺度参数及取值表
尺度参数σ 第一尺度 第二尺度 第三尺度 第四尺度 第五尺度 第六尺度
参数值 20 40 60 80 100 120
MSR算法中尺度参数σ的第二种设计方案：为了弥补第一种方案的不足，设置了第二种方案，参数的个数为9，而且数值设置比较大，具体如表2所示；尺度参数数值较大时，去雾霾效果相对比较好，但是细节信息会丢失多一些；
表2 9个尺度参数及取值表
Figure BDA0002392797210000181
两种方案结果的融合：在实际的操作过程中，为达到既去除雾霾又能很好保持细节信息的目的，通常把两种参数设置方案分别执行处理，然后把各自获得的处理结果进行简单的融合处理，即求和取均值，然后把其作为最终结果。
所述的步骤10中进行局部直方图均衡化处理的具体步骤为：
用MSR算法对图像进行雾霾去除，虽然能够获得比较好的效果，但是处理后整体色调有点偏暗；如果接着用直方图均衡化进行处理，可以弥补这个不足；因为直方图处理图像，能调整图像中灰度的分布状况，使处理后的色调偏亮；因此，在本发明中，把Retinex理论和直方图处理进行级联结合，可实现优势互补，达到最佳效果；
进行的直方图均衡化处理，不是进行的整幅图像的全局均衡处理，而是进行的局部均衡化处理，取得的效果更好；
用直方图原理对图像进行处理时，为避免颜色失真，先把RGB模型图像转换成HSI模型图像，然后对HSI模型中的I分量进行处理；处理结束后，再把HIS模型图像转换成RGB图像；
RGB图像和HSI模型图像转换过程如下：
由RGB图像获取HSI图像的H、S和I分量图的计算公式如下
Figure BDA0002392797210000191
其中θ按式(11)计算，
Figure BDA0002392797210000192
Figure BDA0002392797210000193
Figure BDA0002392797210000194
由HSI图像转换成RBG图像的具体过程如下：
如果0°≤H＜120°，则R、G、B分量由式(14)获取；
Figure BDA0002392797210000201
如果120°≤H＜240°，则R、G、B分量由式(15)进行计算，
Figure BDA0002392797210000202
如果240°≤H＜360°，则R、G、B分量由式(16)进行计算，
Figure BDA0002392797210000203
所述的步骤11中进行小波多尺度分解和融合处理的具体步骤为：
步骤11.1：采用小波多尺度分析理论对前面获得的相位幅频特征图和雾霾去除结果图进行融合处理，得到最后的结果，既实现了雾霾去除，又保持了丰富的细节信息和清晰的对比度。本发明中采用二维平稳小波对图像进行融合处理；
步骤11.2：低频系数的融合，用二维平稳小波分别对相位幅度特征图和滤波结果图进行多尺度分解，在各个分解尺度上再分别进行融合处理，把融合后的结果进行小波逆变换，就获得了最后的结果；融合规则的确定是用小波多尺度分解理论进行特征融合的关键技术。在本发明中，采用像素级进行融合处理策略，两幅图像被分解后，分别获得了若干个高频系数和低频系数子图像；对于低频系数子图像，采取系数加权的形式进行融合处理；对于高频系数子图像采取抓大舍小的原则进行处理；设I1(m,n)和I2(m,n)分别表示两幅按不同方法处理后的遥感图像，IF(m,n)表示融合后的图像，低频系数的加权融合规则数学模型如式(17)所示，
IF(m,n)＝w1·I1(m,n)+w2·I2(m,n) (17)
式中的w1和w2分别表示图像I1(m,n)和图像I2(m,n)的融合权值，即加权系数，而且w1+w2＝1；如果w1＝w2＝0.5，则是平均加权，即均值处理；利用加权系数可以很好地去除源图像中的冗余信息，对于差异性不大的源图像，能够取得很好的效果。在本发明中，两种不同方式处理的特征图来源于同一幅遥感图像，所以他们的差异非常小，用这种方式融合会比其他方式取得更好的效果；
步骤11.3：高频系数的融合，对于每个分解尺度上高频系数子图像的融合采用取大放小的原则进行融合处理，其融合规则数学模型如式(18)所示，
Figure BDA0002392797210000211
这种取极值方法，由于忽略了另外一幅源图像的信息，如果两幅图像差异比较大，融合效果通常不太理想。在本发明中，针对两种方案处理的结果进行融合处理，所以不论是低频系数融合还是高频系数融合，均能达到好的效果。
所述的步骤12中输出结果的具体为：根据获取的图像色彩信息，如果最初输入的是灰度图像，还需转换成灰度图像输出；如果是彩色图像，则直接输出相应的结果。
参照图2，图2中A-J分别代表不同的城市遥感原始雾霾图像，a-e分别代表不同方法获得的实验结果，对比本发明方法和暗通道先验法(DCP)方法、亮度保持模糊直方图(BPDFHE)方法、直方图均衡化(HE)方法、多尺度模型及直方图特性(MSMHC)方法，能够明显看出DCP方法对遥感图像中雾霾的处理效果不是理想，BPDFHE方法能保持较好的亮度，但是去雾霾能力不强，HE方法通常能有效去除遥感图像中的雾霾，其要求图像中雾霾分布比较均匀而且是淡的雾霾，易引起色彩失真，MSMHC算法有效去除雾霾，但在边缘处或灰度变化剧烈的地方容易产生光晕现象，本发明方法不仅能有效去除遥感图像中的雾霾，而且能得到对比度增强的细节信息。
本发明的创新体现在以下几个方面：
(1)根据城市遥感图像的特点：具有丰富的纹理特征，而且纹理特征和地物目标的空间分布不受雾霾的影响，于是提出利用遥感图像的相位一致性特征来增强遥感图像的边缘和细节信息，使遥感图像获得更清晰更丰富的增强细节信息。
(2)相位一致性特征是由相位信息和幅度信息融合获得，不是简单提取相位信息，实现幅度信息和相位信息的卷积及融合。
(3)为了确保遥感图像中雾霾的有效去除，根据多尺度Retinex理论和直方图统计特征各自的优势，设计了两者结合的级联去雾霾策略，做到尽量减少了颜色偏移和细节信息损失。Retinex理论处理图像后，会对图像灰度动态范围进行压缩，使整体颜色偏灰，易出现光晕现象。直方图均衡化处理图像后，会对图像的灰度范围进行拉伸，使整体颜色偏亮。它们两者结合进行级联处理，实现优势互补，尽量去除了遥感图像中的雾霾，还能尽可能保持图像色调。
(4)用Retinex理论处理图像雾霾时，设计了两种不同尺度参数数目和数值的多尺度处理策略，不仅增加了算法的适应范围，而且能确保算法的有效性和特征的逼原性。
(5)采用多尺度平稳小波变换理论完成幅相特征图和雾霾处理后的图像的融合处理，确保融合的精确度，使恢复的信息更准确、更完整。
(6)具有较强的普适性，可以处理灰度和多光谱遥感图像。

Claims (10)

1.一种去除城市遥感图像中雾霾的方法，其特征在于，包括以下步骤：
步骤1：输入城市遥感图像；
步骤2：获取城市遥感图像的色彩信息；
步骤3：把色彩信息中单波段的灰色图像转换成彩色图像；
步骤4：对输入的城市遥感图像进行傅里叶变换；
步骤5：获取相位特征；
步骤6：构造对数Gabor滤波器并提取图像幅度特征；
步骤7：对幅度特征进行颜色补偿；
步骤8：获得相位幅度特征；
步骤9：用多尺度Retinex理论对城市遥感图像进行雾霾去除处理；
步骤10：进行局部直方图均衡化处理；
步骤11：进行小波多尺度分解和融合处理；
步骤12：输出结果。
2.根据权利要求1所述的一种去除城市遥感图像中雾霾的方法，其特征在于：所述的步骤1中输入的遥感图像是单波段灰色图像，或是多光谱图像或多个波段合成的彩色图像，或是由CCD相机获得的RGB真彩色图像，它们都是在雾霾存在的条件下获得的城市区域遥感图像。
3.根据权利要求1所述的一种去除城市遥感图像中雾霾的方法，其特征在于，所述的步骤3中把单波段的灰色图像转换成彩色图像的具体步骤为；
对输入的全色灰度遥感图像进行彩色增强处理，输入的单波段灰度遥感图像用I(m,n)表示，转换后的RGB假彩色遥感图像用IC(m,n)表示，转换的过程如式(1)所示，就是把灰度图像I(m,n)的各个像素值同时分别赋予给假彩色RGB图像IC(m,n)的三个通道，即ICR(m,n)通道，ICG(m,n)通道和ICB(m,n)通道；
Figure FDA0002392797200000021
4.根据权利要求1所述的一种去除城市遥感图像中雾霾的方法，其特征在于，所述的步骤6中构造对数Gabor滤波器并提取图像幅度特征的具体步骤为：
步骤6.1：在频率域直接获取每个像素的幅度值，遥感图像幅度特征的提取与滤波是同时进行，在频率域中，通过构造对数Gabor滤波器对频域中的幅度特征进行滤波，把获得的结果进行傅里叶逆变换，就获得了遥感图像的幅度特征；
步骤6.2：构造Gabor滤波器，对Gabor滤波器设计5个尺度，每个尺度设计6个方向，总共有30个滤波器，把幅度特征经滤波器处理，获得30幅特征图；
步骤6.3：把获得的30幅特征图与步骤5获得的相位信息分别进行卷积运算处理，然后把不同尺度不同方向的处理结果进行加权求和，这里设置各权值相等，即进行求和后再取均值运算，结果作为最终的幅度特征。
5.根据权利要求3所述的一种去除城市遥感图像中雾霾的方法，其特征在于，所述的步骤7中对幅度特征进行颜色补偿的具体步骤为：
利用步骤2获取的颜色信息对幅度特征进行色彩补偿处理，处理的模型公式如下：
Figure FDA0002392797200000031
Figure FDA0002392797200000032
式中FAi(m,n)和
Figure FDA0002392797200000033
分别表示补偿前后的幅度特征，i表示颜色通道数，Ci(m,n)表示第i个颜色通道的色彩补偿系数，β是一个补偿常数；在其他条件不变的情况下，Ci(m,n)值越大，图像颜色失真越大，同时图像偏暗；当其值较小时，颜色失真小，图像偏亮；Ci(m,n)＝20。
6.根据权利要求5所述的一种去除城市遥感图像中雾霾的方法，其特征在于，所述的步骤8中获得相位幅度特征的具体步骤为：
按式(4)把相位特征和幅度特征进行融合处理，以获得最终的相位幅度特征图，
FAP(m,n)＝α·FA(m,n)+β·FP(m,n) (4)
式中FAP(m,n)表示融合后的幅度相位特征，FA(m,n)表示幅度特征，FP(m,n)表示相位特征，α和β表示加权系数，分别设为0.45和0.55。
7.根据权利要求6所述的一种去除城市遥感图像中雾霾的方法，其特征在于，所述的步骤9中用多尺度Retinex理论对城市遥感图像进行雾霾去除处理的具体步骤为：
步骤9.1：获取的图像I(m,n)看作是环境入射光和物体反射相乘的结果，那么它可以用式(5)来表示。
I(m,n)＝L(m,n)·R(m,n) (5)
式中I(m,n)表示用光学传感器获取的图像，即受到了雾霾影响而降质的图像，(m,n)表示像素点的空间位置；L(m,n)表示入射的环境光，即要消除的干扰信息；R(m,n)表示物体表面反射光，即需保留的图像信息；
对公式(5)两边取对数，并进行调整，得到式(6)，
r(m,n)＝log[R(m,n)]＝log[I(m,n)]-log[L(m,n)] (6)
其中L(m,n)用下式进行计算，
L(m,n)＝G(m,n)*I(m,n) (7)
式中的G(m,n)为高斯函数，其定义如下：
Figure FDA0002392797200000041
式中参数σ既为高斯函数的标准偏差，又为其尺度参数；
多尺度Retinex(MSR)算法就是把不同单个尺度Retinex(SSR)算法处理结果进行加权求和，然后把它作为最后的处理结果，其定义如式(9)所示，
Figure FDA0002392797200000042
式中的i为颜色通道数，i∈{1,2,3}，对应于彩色图像的R通道、G通道和B通道；M表示尺度参数值的个数，j表示参数值M中的第j个尺度参数σ，rji(m,n)表示第j尺度下第i颜色通道的反射分量，ri(m,n)表示第i颜色通道经Retinex算法处理后的反射分量，wj为尺度参数值的权值；
步骤9.2：MSR算法中尺度参数σ的个数M及其数值的设计，设置了两种多尺度参数个数和数值方案，具体如下：
MSR算法中尺度参数σ的第一种设计方案：此方案中只设置了6个尺度参数，第一尺度20、第二尺度40、第三尺度60、第四尺度80、第五尺度100、第六尺度120；
MSR算法中尺度参数σ的第二种设计方案：此方案中只设置了9个尺度参数，第一尺度32、第二尺度64、第三尺度128、第四尺度256、第五尺度512、第六尺度1024、第七尺度2048、第八尺度4096、第九尺度8192；
两种方案结果的融合：把两种参数设置方案分别执行处理，然后把各自获得的处理结果进行简单的融合处理，即求和取均值，然后把其作为最终结果。
8.根据权利要求7所述的一种去除城市遥感图像中雾霾的方法，其特征在于，所述的步骤10中进行局部直方图均衡化处理的具体步骤为：
用MSR算法对图像进行雾霾去除，处理后整体色调偏暗；接着用直方图均衡化进行处理；用直方图原理对图像进行处理时，先把RGB模型图像转换成HSI模型图像，然后对HSI模型中的I分量进行处理；处理结束后，再把HIS模型图像转换成RGB图像；
RGB图像和HSI模型图像转换过程如下：
由RGB图像获取HSI图像的H、S和I分量图的计算公式如下
Figure FDA0002392797200000061
其中θ按式(11)计算，
Figure FDA0002392797200000062
Figure FDA0002392797200000063
Figure FDA0002392797200000064
由HSI图像转换成RBG图像的具体过程如下：
如果0°≤H＜120°，则R、G、B分量由式(14)获取；
Figure FDA0002392797200000065
如果120°≤H＜240°，则R、G、B分量由式(15)进行计算，
Figure FDA0002392797200000066
如果240°≤H＜360°，则R、G、B分量由式(16)进行计算，
Figure FDA0002392797200000071
9.根据权利要求8所述的一种去除城市遥感图像中雾霾的方法，其特征在于，所述的步骤11中进行小波多尺度分解和融合处理的具体步骤为：
步骤11.1：采用二维平稳小波对图像进行融合处理；
步骤11.2：低频系数的融合，用二维平稳小波分别对相位幅度特征图和滤波结果图进行多尺度分解，在各个分解尺度上再分别进行融合处理，把融合后的结果进行小波逆变换，就获得了最后的结果；采用像素级进行融合处理策略，两幅图像被分解后，分别获得了若干个高频系数和低频系数子图像；对于低频系数子图像，采取系数加权的形式进行融合处理；对于高频系数子图像采取抓大舍小的原则进行处理；设I1(m,n)和I2(m,n)分别表示两幅按不同方法处理后的遥感图像，IF(m,n)表示融合后的图像，低频系数的加权融合规则数学模型如式(17)所示，
IF(m,n)＝w1·I1(m,n)+w2·I2(m,n) (17)
式中的w1和w2分别表示图像I1(m,n)和图像I2(m,n)的融合权值，即加权系数，而且w1+w2＝1；如果w1＝w2＝0.5，则是平均加权，即均值处理；
步骤11.3：高频系数的融合，对于每个分解尺度上高频系数子图像的融合采用取大放小的原则进行融合处理，其融合规则数学模型如式(18)所示，
Figure FDA0002392797200000081
10.根据权利要求1所述的一种去除城市遥感图像中雾霾的方法，其特征在于，所述的步骤12中输出结果的具体为：根据获取的图像色彩信息，如果最初输入的是灰度图像，还需转换成灰度图像输出；如果是彩色图像，则直接输出相应的结果。
CN202010120431.7A 2020-02-26 2020-02-26 一种去除城市遥感图像中雾霾的方法 Active CN111275652B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010120431.7A CN111275652B (zh) 2020-02-26 2020-02-26 一种去除城市遥感图像中雾霾的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010120431.7A CN111275652B (zh) 2020-02-26 2020-02-26 一种去除城市遥感图像中雾霾的方法

Publications (2)

Publication Number Publication Date
CN111275652A true CN111275652A (zh) 2020-06-12
CN111275652B CN111275652B (zh) 2023-03-10

Family

ID=71004069

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010120431.7A Active CN111275652B (zh) 2020-02-26 2020-02-26 一种去除城市遥感图像中雾霾的方法

Country Status (1)

Country Link
CN (1) CN111275652B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112001288A (zh) * 2020-08-17 2020-11-27 西京学院 一种单幅多光谱遥感图像深灰色飞机快速检测方法
CN112785521A (zh) * 2021-01-19 2021-05-11 澜途集思生态科技集团有限公司 一种雾霾条件下的遥感影像处理方法
CN112907460A (zh) * 2021-01-25 2021-06-04 宁波市鄞州区测绘院 一种遥感图像增强方法
CN113160073A (zh) * 2021-03-19 2021-07-23 西安邮电大学 一种滚动深度学习与Retinex理论结合的遥感图像雾霾去除方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9465994B1 (en) * 2015-02-23 2016-10-11 Amazon Technologies, Inc. Predicting performance and success of large-scale vision algorithms
CN109003238A (zh) * 2018-06-26 2018-12-14 西京学院 一种基于模型与直方图及灰度增强的图像雾霾去除方法
CN109447917A (zh) * 2018-10-31 2019-03-08 西京学院 基于内容和特征及多尺度模型的遥感图像雾霾消除方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9465994B1 (en) * 2015-02-23 2016-10-11 Amazon Technologies, Inc. Predicting performance and success of large-scale vision algorithms
CN109003238A (zh) * 2018-06-26 2018-12-14 西京学院 一种基于模型与直方图及灰度增强的图像雾霾去除方法
CN109447917A (zh) * 2018-10-31 2019-03-08 西京学院 基于内容和特征及多尺度模型的遥感图像雾霾消除方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
王荣荣等: "基于复Shearlet变换耦合改进引导滤波的图像融合算法", 《电子测量与仪器学报》 *
胡根生等: "基于对偶树复小波变换的遥感图像薄云去除", 《安徽大学学报(自然科学版)》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112001288A (zh) * 2020-08-17 2020-11-27 西京学院 一种单幅多光谱遥感图像深灰色飞机快速检测方法
CN112001288B (zh) * 2020-08-17 2024-02-20 广州商学院 一种单幅多光谱遥感图像深灰色飞机快速检测方法
CN112785521A (zh) * 2021-01-19 2021-05-11 澜途集思生态科技集团有限公司 一种雾霾条件下的遥感影像处理方法
CN112907460A (zh) * 2021-01-25 2021-06-04 宁波市鄞州区测绘院 一种遥感图像增强方法
CN112907460B (zh) * 2021-01-25 2022-07-29 宁波市鄞州区测绘院 一种遥感图像增强方法
CN113160073A (zh) * 2021-03-19 2021-07-23 西安邮电大学 一种滚动深度学习与Retinex理论结合的遥感图像雾霾去除方法
CN113160073B (zh) * 2021-03-19 2023-10-03 西安邮电大学 一种滚动深度学习与Retinex理论结合的遥感图像雾霾去除方法

Also Published As

Publication number Publication date
CN111275652B (zh) 2023-03-10

Similar Documents

Publication Publication Date Title
CN111275652B (zh) 一种去除城市遥感图像中雾霾的方法
CN112233026A (zh) 一种基于多尺度残差注意力网络的sar图像去噪方法
CN109447917B (zh) 基于内容和特征及多尺度模型的遥感图像雾霾消除方法
CN110400267A (zh) 一种基于巡检图像的预处理方法
Guo et al. Automatic image haze removal based on luminance component
CN108765336A (zh) 基于暗亮原色先验与自适应参数优化的图像去雾方法
CN112070717B (zh) 基于图像处理的输电线路覆冰厚度检测方法
CN107784639B (zh) 一种无人机遥感图像改进多边滤波去噪方法
CN113313639A (zh) 基于Retinex多层次分解的图像增强方法
CN115578660B (zh) 基于遥感图像的土地地块分割方法
CN109003238B (zh) 一种基于模型与直方图及灰度增强的图像雾霾去除方法
CN112541869A (zh) 一种基于matlab实现的Retinex图像去雾方法
Wen et al. Autonomous robot navigation using Retinex algorithm for multiscale image adaptability in low-light environment
Zhang et al. Salient target detection based on the combination of super-pixel and statistical saliency feature analysis for remote sensing images
CN117575953B (zh) 一种高分辨率林业遥感图像细节增强方法
Fu et al. Scene-awareness based single image dehazing technique via automatic estimation of sky area
CN113822816A (zh) 气雾散射模型优化的单张遥感图像去霾方法
CN111612797B (zh) 一种水稻图像信息处理***
CN113129300A (zh) 一种降低误检率的排水管道缺陷检测方法、装置、设备及介质
CN112488961A (zh) 一种基于对数方程的t参数同态滤波方法
CN115393406B (zh) 一种基于孪生卷积网络的图像配准方法
CN111461999A (zh) 一种基于超像素相似性测量的sar图像相干斑抑制方法
CN115908155A (zh) Nsst域结合gan及尺度相关系数的低照度图像增强及去噪方法
Khalaf et al. Hazy Image Enhancement Using DCP and AHE Algorithms with YIQ Color Space
Zheng An exploration of color fusion with multispectral images for night vision enhancement

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