CN113570536B - 基于cpu和gpu协同处理的全色和多光谱影像实时融合方法 - Google Patents

基于cpu和gpu协同处理的全色和多光谱影像实时融合方法 Download PDF

Info

Publication number
CN113570536B
CN113570536B CN202110877331.3A CN202110877331A CN113570536B CN 113570536 B CN113570536 B CN 113570536B CN 202110877331 A CN202110877331 A CN 202110877331A CN 113570536 B CN113570536 B CN 113570536B
Authority
CN
China
Prior art keywords
image
panchromatic
multispectral
virtual
fusion
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
CN202110877331.3A
Other languages
English (en)
Other versions
CN113570536A (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.)
No61646 Unit Of Pla
Original Assignee
No61646 Unit Of Pla
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 No61646 Unit Of Pla filed Critical No61646 Unit Of Pla
Priority to CN202110877331.3A priority Critical patent/CN113570536B/zh
Publication of CN113570536A publication Critical patent/CN113570536A/zh
Application granted granted Critical
Publication of CN113570536B publication Critical patent/CN113570536B/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
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
    • 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
    • G06T2207/10036—Multispectral image; Hyperspectral 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
    • G06T2207/10041—Panchromatic 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/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging

Landscapes

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

Abstract

本发明公开了一种基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，对全色多光谱数据进行快速融合处理，生成具有高空间分辨率及高光谱分辨率的融合影像。本发明将全色与多光谱融合的所有步骤都在内存中完成，并且将计算量较大的步骤都映射到GPU进行，能够在保障融合精度的同时极大地提高融合效率，通过RFM考虑了全色多光谱之间严密的几何对应关系，再通过微分纠正最大限度地限制了全色多光谱之间可能存在的几何偏差，使得全色与多光谱的配准效果能够达到最优，最后通过基于多尺度的SFIM融合方法与全色光谱分解的影像融合方法获取效果最优的融合图像。本发明适用于大数据量的传感器数据快速融合处理需求。

Description

基于CPU和GPU协同处理的全色和多光谱影像实时融合方法
技术领域
本发明属于图像处理领域，具体涉及一种基于CPU和GPU协同处理的全色和多光谱影像实时融合方法。
背景技术
从传感器的工作方式的角度看，影像数据分为多光谱影像和全色波段图像。全色波段图像数据具有比较高的空间分辨率，而光谱分辨率较低；多光谱影像数据由于采用窄波段光源进行成像，光谱分辨率较高而空间分辨率很低。全色多光谱融合方法就是将这两者的优点结合起来，使融合后图像既具有高的空间分辨能力又具有高的光谱分辨能力，以利于后期的判读和处理。
对于相机而言，探测像元，简称探元，只有接受到足够的能量才能维持信噪比，采用多光谱单波段成像时，传感器接收能量弱，只能增大探元尺寸，从而导致多光谱影像的空间分辨率低于全色波段图像；而全色波段图像的空间分辨率虽然高于多光谱，但是为了增加接收能量的强度，其传感器波谱范围较宽，其图像的光谱分辨能力低于多光谱。在相等的条件下，提高空间分辨率必然会降低光谱分辨率，这是不可调和的矛盾。因此在相机实际设计的过程中，高频空间信息与光谱信息相分离，一般的光学传感器只提供高分辨率全色影像和低分辨率多光谱影像。如何得到光谱信息与多光谱影像的融合影像是研究的重点。
随着社会、经济和科技的发展，影像分辨率已经达到米级或亚米级，其单景一级全色影像数据量也逐步增加，全色数据量的增加也会导致融合数据的迅速增加，最高可达上百GB。当任务对时效性的要求非常高，甚至到分钟级响应需求时，对光谱图像的融合处理算法提出了很高的要求，传统的融合处理方法延时较高，难以满足快速处理的需求。
发明内容
为了突破全色多光谱快速融合的瓶颈，满足不同任务快速反应的重大需求，迫切需要研究新的全色多光谱影像融合高性能处理的方法，本发明公开了一种基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，在保障融合效果的前提下，有效地提升了大数据量融合的效率，具有一定理论价值与实际意义。
为了解决大数据量的全色和多光谱融合过程中效率较低的问题，本发明主要基于CPU/GPU协同处理，整个处理流程都在内存当中进行，CPU负责流程控制与并行规则划分，GPU负责并行处理。本发明提出了一种基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，其具体步骤包括：
S101，通过有理函数模型计算全色影像和多光谱影像的重叠区域；
通过有理函数模型RFM确定全色影像与多光谱影像各自的物方坐标，物方坐标表达式为(B,L,H)，其包括大地纬度B、大地经度L、大地高H，H为平均的大地高程。
对物方坐标进行归一化处理，得到归一化后的物方坐标(U,V,W)，其计算公式为：
Figure BDA0003190811590000021
其中，LonOff，LatOff，HeiOff分别为物方坐标(B,L,H)的平移值，LonScale，LatScale，HeiScale分别为物方坐标(B,L,H)的缩放值。
对于每一景影像，归一化后的物方坐标(U,V,W)与像方坐标(sn,ln)的关系用多项式比值形式表示为：
Figure BDA0003190811590000022
其中，等号右边的多项式分子、分母的表达式分别为：
NumL(U，V，W)＝a1+a2V+a3U+a4W+a5VU+a6VW+a7UW+a8V2+a9U2+a10W2+a11VUW+a12V3+a13VU2+a14VW2+a15V2U+a16U3+a17UW2+a18V2W+a19U2W+a20W3，
DenL(U，V，W)＝b1+b2V+b3U+b4W+b5VU+b6VW+b7UW+b8V2+b9U2+b10W2+b11VUW+b12V3+b13VU2+b14VW2+b15V2U+b16U3+b17UW2+b18V2W+b19U2W+b20W3，
NumS(U，V，W)＝c1+c2V+c3U+c4W+c5VU+c6VW+c7UW+c8V2+c9U2+c10W2+c11VUW+c12V3+c13VU2+c14VW2+c15V2U+c16U3+c17UW2+c18V2W+c19U2W+c20W3，
DenS(U，V，W)＝d1+d2V+d3U+d4W+d5VU+d6VW+d7UW+d8V2+d9U2+d10W2+d11VUW+d12V3+d13VU2+d14VW2+d15V2U+d16U3+d17UW2+d18V2W+d19U2W+d20W3，
其中，ai，bi，ci，di，i＝1，2，…，20，为有理多项式系数(Rational PolynomialCoefficients，RPCs)。
得到全色影像与多光谱影像的地理范围之后，计算出全色影像与多光谱影像的重叠区域的地理范围，并且分别将两张影像各自重叠区域的地理范围通过RFM模型，得到其在各自影像的像方坐标系下的坐标值。
S102，读取全色影像与多光谱影像的重叠区域中的全色影像，并且以其为基准，通过GPU对多光谱影像进行虚拟化，得到虚拟化后的多光谱影像，即虚拟影像；
将步骤S101中计算得到的全色影像与多光谱影像的重叠区域中的全色影像读入CPU内存当中，以该重叠区域中的全色影像为基准，通过GPU对多光谱影像进行虚拟化，其具体步骤包括：
S1021，以全色影像的像方坐标为基准，将虚拟影像划分为均匀的规则格网，并且计算每个格网的四个角点的像方坐标；
S1022，对虚拟影像的每个格网的四个角点的像方坐标(s，l)，通过RFM模型计算其对应的物方坐标；
S1023，对虚拟影像的每个格网的四个角点的物方坐标，通过RFM模型计算多光谱影像中对应点的像方坐标(s′，l′)；
S1024，利用虚拟影像每个格网的四个角点的像方坐标，构建虚拟影像像方坐标与多光谱影像像方坐标的变换模型g()，其变换过程表示为(s，l)＝g(s′，l′)，该变换模型g()采用线性变换模型、仿射变换模型或透视变换模型中的一种。
S1025，对虚拟影像像素进行重采样；将虚拟影像像素重采样的过程映射到GPU进行处理，将GPU计算过程中的每一个CUDA线程设置为仅对一个像素进行计算。具体来说，将虚拟影像上的每个像素点的坐标通过变换模型g()，计算得到多光谱影像上的像素点坐标，并且选用最近邻近采样、双线性采样或双三次重采样方法中的一种。
S103，将虚拟化后的多光谱影像转成虚拟全色影像，通过GPU将虚拟全色影像与真实的全色影像进行格网点匹配，得到同名网格点；
首先将虚拟化后的多光谱影像中的每个像素转换成与全色影像光学特性一致的虚拟全色单一波段像素，该过程表示为：
Mssp＝e1*BandB+e2*BandG+e3*BandR+e4*BandNIR，
上式中，Mssp为转换后的与全色影像光学特性一致的虚拟全色单一波段像素值，e1、e2、e3、e4分别为蓝光波段B、绿光波段G、红光波段R、近红外波段NIR对应的波段范围占全色波段范围的比例，BandB、BandG、BandR、BandNIR分别代表多光谱虚拟影像中的蓝光波段、绿光波段、红光波段以及近红外波段影像的像素值。
提取全色影像与虚拟全色影像的同名点，该同名点作为多光谱虚拟影像进行数字微分纠正的输入。所述的提取全色影像与虚拟全色影像的同名点，采用相关系数法，以全色影像与虚拟全色影像的相关系数作为影像匹配的匹配测度，在指定的影像搜索范围中计算相关系数，而后选取相关系数最大的点作为同名点，相关系数是标准化的协方差函数。
选择M波段的全色影像上的某一像素点为初值点，以该初值点为中心，建立m×n个像素的匹配窗口，并在N波段的虚拟全色影像上建立与M波段的全色影像的初值点上同样大小和形状的匹配窗口，计算两个匹配窗口的相关系数，其计算公式为：
Figure BDA0003190811590000051
其中，R(c,r)表示两个匹配窗口在行和列上的坐标差分别为c和r时的相关系数，m,n分别表示匹配窗口的行和列大小，i,j分别表示匹配窗口内的行序号和列序号，gi,j表示M波段的全色影像和N波段的虚拟全色影像的匹配窗口内第i行、第j列的灰度值，gi+r,j+c表示上述的两个影像的匹配窗口内第i+r行、第j+c列的灰度值，c,r分别是N波段的虚拟全色影像的匹配窗口对于M波段的全色影像的匹配窗口中的同名点在行、列上的坐标差。对于已经进行了格网点匹配处理的多光谱影像，取c＝r＝0；对于未进行格网点匹配的影像，c和r的取值则根据相机的硬件参数值来确定。
在全色影像上从设定的格网点中确定控制点，采用相关系数法，依次获取每个格网点在虚拟全色影像上的同名点。将相关系数的计算过程映射到GPU上进行，GPU上的每个线程块对应一个格网点，线程块中的每个线程实现一组相关系数的计算，每个线程块计算完毕之后获取相关系数最大的同名点，并且利用抛物线法对相关系数最大的同名点相邻的四个点进行拟合，得到该相关系数最大的同名点的亚像素值。
S104，获取同名格网点之后，通过GPU对多光谱虚拟影像进行数字微分纠正，得到数字微分纠正后的多光谱影像；
对多光谱虚拟影像进行数字微分纠正，多光谱虚拟影像中的每个格网包括四组同名点，数字微分纠正过程f(x'，y')的数学表达式为：
Figure BDA0003190811590000052
式中，(x,y)为所述的同名点在数字微分纠正前的像方坐标，(x',y')为所述的同名点在数字微分纠正后的像方坐标，g1、f1、g2、f2、g3、f3、g4、f4分别为数字微分纠正过程的8个变换系数，该8个变换系数通过对4组同名点联合建立的8个数字微分纠正过程方程进行解算得到。分别解算每个格网的变换系数，从而对多光谱虚拟影像进行数字微分纠正，再进行最近邻近采样、双线性采样或双三次重采样。将数字微分纠正的过程映射到GPU进行处理，将GPU计算过程中的每一个CUDA线程设置为仅对一个像素进行计算。
S105，对数字微分纠正后的多光谱影像以及全色影像进行快速融合，得到最终的融合结果。
所述的步骤S105，通过两种方式实施，一是利用基于多尺度的平滑滤波亮度变换模型SFIM的融合算法(MSSFIM)实现，二是利用基于全色光谱分解的影像融合方法(PSD)实现。
所述的基于多尺度的SFIM模型的融合算法，是基于多尺度原理与SFIM模型的影像融合算法；SFIM模型的融合过程表示为，
Figure BDA0003190811590000061
其中，MS是多光谱影像的像素值，Pan是全色影像的像素值，Pan′是低分辨率全色影像像素值，Fusion是融合结果。Pan和Pan′的比值只保留了高分辨率图像的边缘细节信息。
对全色影像利用高斯金字塔方法，生成不同尺度的低分辨率全色影像，并向SFIM模型得到的融合影像中添加不同清晰程度的高频信息，以增强融合后的图像的清晰度。
所述的基于全色光谱分解的影像融合方法，利用了同时拍摄的全色影像与多光谱影像之间的强相关性，通过线性变换加残差消减的方式将全色影像像元亮度值DN分解为多光谱影像波段像元亮度值DN，再利用比值变换法进行影像融合。在比值变换法中，全色影像与所生成不同尺度的低分辨率全色影像的像素值比值等于融合图像与低质量多光谱影像的像素值的比值，即：
Figure BDA0003190811590000071
其中，PanHQ为高质量全色影像，即原始全色影像，PanHQ(i，j)表示原始全色影像的第i行、第j列像素点的像素值；PanLQ为所生成不同尺度的低质量全色影像，PanLQ(i，j)表示所生成不同尺度的低质量全色影像的第i行、第j列像素点的像素值；
Figure BDA0003190811590000072
为理想的高质量多光谱影像，即融合图像，
Figure BDA0003190811590000073
表示融合图像的g波段图像第i行、第j列像素点的像素值；
Figure BDA0003190811590000074
为低质量多光谱影像，其通过对原始多光谱影像进行采样得到，
Figure BDA0003190811590000075
表示低质量多光谱影像的g波段图像第i行、第j列像素点的像素值；i，j，g分别是行、列、波段的序号。全色影像与多光谱影像之间的数学关系描述为：
Pan＝k(g)MS(g)+b(g)T+E(g)，
其中，b(g)是多光谱影像第g波段对应的偏置系数，k(g)为多光谱影像第g波段对应的比例系数，T为维度与Pan相同、且其元素都为1的矩阵，E(g)是多光谱影像第g波段对应的残差矩阵，Pan是全色影像的像素值，MS(g)是g波段多光谱影像的像素值。通过低分辨率的全色多光谱影像拟合得到比例系数、偏置系数及残差矩阵，并利用这些系数得到融合的高分辨率全色影像。
本发明的有益效果为：
相对于传统的分步骤融合方法，本方法将全色与多光谱融合的所有步骤都在内存中完成，并且将计算量较大的步骤都映射到GPU进行，能够在保障融合精度的同时极大地提高融合效率。并且该方法不仅通过RFM考虑了全色多光谱之间严密的几何对应关系，而且还通过微分纠正最大限度地限制了全色多光谱之间可能存在的几何偏差，使得全色与多光谱的配准效果能够达到最优。最后该方法通过基于多尺度的SFIM融合方法与全色光谱分解的影像融合方法获取效果最优的融合图像。本发明适用于大数据量的传感器影像快速融合处理需求。
附图说明
图1是本发明一种基于CPU/GPU协同处理的全色多光谱近实时融合方法的流程示意图。
图2为全色多光谱的原始影像和虚拟影像坐标映射流程图。
具体实施方式
为了更好的了解本发明内容，这里给出一个实施例。
本发明公开了一种基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，在保障融合效果的前提下，有效地提升了大数据量融合的效率，具有一定理论价值与实际意义。图1是本发明一种基于CPU/GPU协同处理的全色多光谱近实时融合方法的流程示意图。
为了解决大数据量的全色和多光谱融合过程中效率较低的问题，本发明主要基于CPU/GPU协同处理，整个处理流程都在内存当中进行，CPU负责流程控制与并行规则划分，GPU负责并行处理。本发明提出了一种基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，其具体步骤包括：
S101，通过有理函数模型计算全色影像和多光谱影像的重叠区域；
通过有理函数模型RFM确定全色影像与多光谱影像各自的物方坐标，物方坐标表达式为(B,L,H)，其包括大地纬度B、大地经度L、大地高H，H为平均的大地高程。
对物方坐标进行归一化处理，得到归一化后的物方坐标(U,V,W)，其计算公式为：
Figure BDA0003190811590000091
其中，LonOff，LatOff，HeiOff分别为物方坐标(B，L，H)的平移值，LonScale，LatScale，HeiScale分别为物方坐标(B，L，H)的缩放值。
对于每一景影像，归一化后的物方坐标(U，V，W)与像方坐标(sn，ln)的关系用多项式比值形式表示为：
Figure BDA0003190811590000092
其中，等号右边的多项式分子、分母的表达式分别为：
NumL(U，V，W)＝a1+a2V+a3U+a4W+a5VU+a6VW+a7UW+a8V2+a9U2+a10W2+a11VUW+a12V3+a13VU2+a14VW2+a15V2U+a16U3+a17UW2+a18V2W+a19U2W+a20W3，
DenL(U，V，W)＝b1+b2V+b3U+b4W+b5VU+b6VW+b7UW+b8V2+b9U2+b10W2+b11VUW+b12V3+b13VU2+b14VW2+b15V2U+b16U3+b17UW2+b18V2W+b19U2W+b20W3，
NumS(U，V，W)＝c1+c2V+c3U+c4W+c5VU+c6VW+c7UW+c8V2+c9U2+c10W2+c11VUW+c12V3+c13VU2+c14VW2+c15V2U+c16U3+c17UW2+c18V2W+c19U2W+c20W3，
DenS(U，V，W)＝d1+d2V+d3U+d4W+d5VU+d6VW+d7UW+d8V2+d9U2+d10W2+d11VUW+d12V3+d13VU2+d14VW2+d15V2U+d16U3+d17UW2+d18V2W+d19U2W+d20W3，
其中，ai，bi，ci，di，i＝1，2，…，20，为有理多项式系数(Rational PolynomialCoefficients，RPCs)。一般情况下，b1，d1均取值为1。
得到全色影像与多光谱影像的地理范围之后，计算出全色影像与多光谱影像的重叠区域的地理范围，并且分别将两张影像各自重叠区域的地理范围通过RFM模型，得到其在各自影像的像方坐标系下的坐标值。
S102，读取全色影像与多光谱影像的重叠区域中的全色影像，并且以其为基准，通过GPU对多光谱影像进行虚拟化，得到虚拟化后的多光谱影像，即虚拟影像；图2为全色多光谱的原始影像和虚拟影像坐标映射流程图。
将步骤S101中计算得到的全色影像与多光谱影像的重叠区域中的全色影像读入CPU内存当中，以该重叠区域中的全色影像为基准，通过GPU对多光谱影像进行虚拟化，其具体步骤包括：
S1021，以全色影像的像方坐标为基准，将虚拟影像划分为均匀的规则格网，并且计算每个格网的四个角点的像方坐标；
S1022，对虚拟影像的每个格网的四个角点的像方坐标(s,l)，通过RFM模型计算其对应的物方坐标；
S1023，对虚拟影像的每个格网的四个角点的物方坐标，通过RFM模型计算多光谱影像中对应点的像方坐标(s'，l')；
S1024，利用虚拟影像每个格网的四个角点的像方坐标，构建虚拟影像像方坐标与多光谱影像像方坐标的变换模型g()，其变换过程表示为(s，l)＝g(s'，l')，该变换模型g()采用线性变换模型、仿射变换模型或透视变换模型中的一种，本发明从精度和效率的角度，优先考虑选用透视变换模型。
S1025，对虚拟影像像素进行重采样；前面所述的步骤已经完成了原始影像与虚拟影像坐标的映射工作，后续将要进行像素的重采样工作，值的一提的是像素重采样是天然并行的工作，为了提高效率，加快处理速度，将虚拟影像像素重采样的过程映射到GPU进行处理，将GPU计算过程中的每一个CUDA线程设置为仅对一个像素进行计算。具体来说，将虚拟影像上的每个像素点的坐标通过变换模型g()，计算得到多光谱影像上的像素点坐标，并且选用最近邻近采样、双线性采样或双三次重采样方法中的一种，本发明从精度和效率的角度考虑优先选用双线性重采样，每次同时处理的像素数即可由不同的GPU的性能所决定。
S103，将虚拟化后的多光谱影像转成虚拟全色影像，通过GPU将虚拟全色影像与真实的全色影像进行格网点匹配，得到同名网格点；
通过步骤S102计算即可在内存中得到全色影像与多光谱的虚拟影像，在步骤S103中，为了更好地与全色影像计算同名点，首先将虚拟化后的多光谱影像中的每个像素转换成与全色影像光学特性一致的虚拟全色单一波段像素，该过程表示为：
Mssp＝e1*BandB+e2*BandG+e3*BandR+e4*BandNIR，
上式中，Mssp为转换后的与全色影像光学特性一致的虚拟全色单一波段像素值，e1、e2、e3、e4分别为蓝光波段B、绿光波段G、红光波段R、近红外波段NIR对应的波段范围占全色波段范围的比例，该比例可由传感器的设计参数得到，BandB、BandG、BandR、BandNIR分别代表多光谱虚拟影像中的蓝光波段、绿光波段、红光波段以及近红外波段影像的像素值。
提取全色影像与虚拟全色影像的同名点，该同名点作为多光谱虚拟影像进行数字微分纠正的输入。所述的提取全色影像与虚拟全色影像的同名点，采用相关系数法，以全色影像与虚拟全色影像的相关系数作为影像匹配的匹配测度，在指定的影像搜索范围中计算相关系数，而后选取相关系数最大的点作为同名点，相关系数是标准化的协方差函数，可以克服灰度的线性变形，保持相关系数的不变性，是影像匹配中最常使用的匹配测度。
选择M波段的全色影像上的某一像素点为初值点，以该初值点为中心，建立m×n个像素的匹配窗口，并在N波段的虚拟全色影像上建立与M波段的全色影像的初值点上同样大小和形状的匹配窗口，计算两个匹配窗口的相关系数，其计算公式为：
Figure BDA0003190811590000121
其中，R(c,r)表示两个匹配窗口在行和列上的坐标差分别为c和r时的相关系数，m,n分别表示匹配窗口的行和列大小，i,j分别表示匹配窗口内的行序号和列序号，gi,j表示M波段的全色影像和N波段的虚拟全色影像的匹配窗口内第i行、第j列的灰度值，gi+r,j+c表示上述的两个影像的匹配窗口内第i+r行、第j+c列的灰度值，c,r分别是N波段的虚拟全色影像的匹配窗口对于M波段的全色影像的匹配窗口中的同名点在行、列上的坐标差。对于已经进行了格网点匹配处理的多光谱影像，取c＝r＝0；对于未进行格网点匹配的影像，c和r的取值则根据相机的硬件参数值来确定。
在全色影像上从设定的格网点中确定控制点，采用相关系数法，依次获取每个格网点在虚拟全色影像上的同名点。每个格网点的相关系数的计算可以并行进行，将相关系数的计算过程映射到GPU上进行，GPU上的每个线程块对应一个格网点，线程块中的每个线程实现一组相关系数的计算，每个线程块计算完毕之后获取相关系数最大的同名点，并且利用抛物线法对相关系数最大的同名点相邻的四个点进行拟合，得到该相关系数最大的同名点的亚像素值。
S104，获取同名格网点之后，通过GPU对多光谱虚拟影像进行数字微分纠正，得到数字微分纠正后的多光谱影像；
对多光谱虚拟影像进行数字微分纠正，多光谱虚拟影像中的每个格网包括四组同名点，数字微分纠正过程f(x'，y')的数学表达式为：
Figure BDA0003190811590000122
式中，(x,y)为所述的同名点在数字微分纠正前的像方坐标，(x',y')为所述的同名点在数字微分纠正后的像方坐标，g1、f1、g2、f2、g3、f3、g4、f4分别为数字微分纠正过程的8个变换系数，该8个变换系数通过对4组同名点联合建立的8个数字微分纠正过程方程进行解算得到。分别解算每个格网的变换系数，从而对多光谱虚拟影像进行数字微分纠正，再进行最近邻近采样、双线性采样或双三次重采样。将数字微分纠正的过程映射到GPU进行处理，将GPU计算过程中的每一个CUDA线程设置为仅对一个像素进行计算。具体来说，将多光谱虚拟影像上的每个点坐标(x，y)通过变换模型f(x′，y′)得到多光谱上的像素坐标(x′，y′)，并且，每次同时处理的效率与像素数即可由不同的GPU的性能所决定。
S105，对数字微分纠正后的多光谱影像以及全色影像进行快速融合，得到最终的融合结果。
所述的步骤S105，通过两种方式实施，一是利用基于多尺度的SFIM模型的融合算法(MSSFIM)实现，二是利用基于全色光谱分解的影像融合方法(PSD)实现。
所述的基于多尺度的SFIM模型的融合算法，是基于多尺度原理与SFIM模型的影像融合算法；改善了原始SFIM模型中图像清晰度不佳的劣势，同时保留了其色彩偏差小的优点。SFIM模型的融合过程表示为，
Figure BDA0003190811590000131
其中，MS是多光谱影像的像素值，Pan是全色影像的像素值，Pan′是低分辨率全色影像像素值，Fusion是融合结果。Pan和Pan′的比值只保留了高分辨率图像的边缘细节信息，而基本消除了光谱与对比度信息。SFIM模型可以被视为是在低分辨率图像中添加了高分辨率细节，其融合结果与高分辨率图像的光谱差异无关。
但是，SFIM模型虽然有较好的色彩保持能力，但是清晰度不佳，图像呈现发虚状态。为此，对全色影像利用高斯金字塔方法，生成不同尺度的低分辨率全色影像，并向SFIM模型得到的融合影像中添加不同清晰程度的高频信息，以增强融合后的图像的清晰度。
所述的基于全色光谱分解的影像融合方法，利用了同时拍摄的全色影像与多光谱影像之间的强相关性，通过线性变换加残差消减的方式将全色影像像元亮度值DN分解为多光谱影像波段像元亮度值DN，再利用比值变换法进行影像融合。同时相的全色与多光谱影像之间存在很强的相关性，在比值变换法中，全色影像与所生成不同尺度的低分辨率全色影像的像素值比值等于融合图像与低质量多光谱影像的像素值的比值，即：
Figure BDA0003190811590000141
其中，PanHQ为高质量全色影像，即原始全色影像，PanHQ(i，j)表示原始全色影像的第i行、第j列像素点的像素值；PanLQ为所生成不同尺度的低质量全色影像，PanLQ(i，j)表示所生成不同尺度的低质量全色影像的第i行、第j列像素点的像素值；
Figure BDA0003190811590000142
为理想的高质量多光谱影像，即融合图像，
Figure BDA0003190811590000143
表示融合图像的g波段图像第i行、第j列像素点的像素值；
Figure BDA0003190811590000144
为低质量多光谱影像，其通过对原始多光谱影像进行采样得到，
Figure BDA0003190811590000145
表示低质量多光谱影像的g波段图像第i行、第j列像素点的像素值；i，j，g分别是行、列、波段的序号。比值变换法将相同的高频信息重复地添加入不同的波段，这是造成融合图像光谱畸变的重要原因。上式成立的基础是全色影像的每个像素与其对应质量的多光谱影像像素都具有无偏置的线性关系，即：
Pan(i，j)＝k(g)*MS(g)(i，j)，
其中，k(g)为第g波段对应的比例系数。该前提在实际中并非完全符合，全色影像与多光谱影像具有一定的发散性。即，全色影像与多光谱影像之间的数学关系描述为：
Pan＝k(g)MS(g)+b(g)T+Ε(g)，
其中，b(g)是多光谱影像第g波段对应的偏置系数，k(g)为多光谱影像第g波段对应的比例系数，T为维度与Pan相同、且其元素都为1的矩阵，E(g)是多光谱影像第g波段对应的残差矩阵。通过低分辨率的全色多光谱影像拟合得到比例系数、偏置系数及残差矩阵，并利用这些系数得到融合的高分辨率全色影像，可以获得效果更好的融合影像。
以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。

Claims (5)

1.一种基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，其特征在于，其具体步骤包括：
S101，通过有理函数模型计算全色影像和多光谱影像的重叠区域；
通过有理函数模型RFM确定全色影像与多光谱影像各自的物方坐标，物方坐标表达式为(B，L，H)，其包括大地纬度B、大地经度L、大地高H，H为平均的大地高程；
对物方坐标进行归一化处理，得到归一化后的物方坐标(U，V，W)，其计算公式为：
Figure FDA0003412361170000011
其中，LonOff，LatOff，HeiOff分别为物方坐标(B，L，H)的平移值，LonScale，LatScale，HeiScale分别为物方坐标(B，L，H)的缩放值；
对于每一景影像，归一化后的物方坐标(U，V，W)与像方坐标(sn，ln)的关系用多项式比值形式表示为：
Figure FDA0003412361170000012
其中，等号右边的多项式分子、分母的表达式分别为：
NumL(U，V，W)＝a1+a2V+a3U+a4W+a5VU+a6VW+a7UW+a8V2+a9U2+a10W2+a11VUW+a12V3+a13VU2+a14VW2+a15V2U+a16U3+a17UW2+a18V2W+a19U2W+a20W3，
DenL(U，V，W)＝b1+b2V+b3U+b4W+b5VU+b6VW+b7UW+b8V2+b9U2+b10W2+b11VUW+b12V3+b13VU2+b14VW2+b15V2U+b16U3+b17UW2+b18V2W+b19U2W+b20W3，
Nums(U，V，W)＝c1+c2V+c3U+c4W+c5VU+c6VW+c7UW+c8V2+c9U2+c10W2+c11VUW+c12V3+c13VU2+c14VW2+c15V2U+c16U3+c17UW2+c18V2W+c19U2W+c20W3，
DenS(U，V，W)＝d1+d2V+d3U+d4W+d5VU+d6VW+d7UW+d8V2+d9U2+d10W2+d11VUW+d12V3+d13VU2+d14VW2+d15V2U+d16U3+d17UW2+d18V2W+d19U2W+d20W3，
其中，ai，bi，ci，di，i＝1，2，…，20，为有理多项式系数；
得到全色影像与多光谱影像的地理范围之后，计算出全色影像与多光谱影像的重叠区域的地理范围，并且分别将两张影像各自重叠区域的地理范围通过RFM模型，得到其在各自影像的像方坐标系下的坐标值；
S102，读取全色影像与多光谱影像的重叠区域中的全色影像，并且以其为基准，通过GPU对多光谱影像进行虚拟化，得到虚拟化后的多光谱影像，即虚拟影像；
S103，将虚拟化后的多光谱影像转成虚拟全色影像，通过GPU将虚拟全色影像与真实的全色影像进行格网点匹配，得到同名网格点；
S104，获取同名格网点之后，通过GPU对多光谱虚拟影像进行数字微分纠正，得到数字微分纠正后的多光谱影像；
S105，对数字微分纠正后的多光谱影像以及全色影像进行快速融合，得到最终的融合结果；
步骤S105利用基于全色光谱分解的影像融合方法实现；
所述的基于全色光谱分解的影像融合方法，利用了同时拍摄的全色影像与多光谱影像之间的强相关性，通过线性变换加残差消减的方式将全色影像像元亮度值DN分解为多光谱影像波段像元亮度值DN，再利用比值变换法进行影像融合；在比值变换法中，全色影像与所生成不同尺度的低分辨率全色影像的像素值比值等于融合图像与低质量多光谱影像的像素值的比值，即：
Figure FDA0003412361170000021
其中，PanHQ为高质量全色影像，即原始全色影像，PanHQ(i，j)表示原始全色影像的第i行、第j列像素点的像素值；PanLQ为所生成不同尺度的低质量全色影像，PanLQ(i，j)表示所生成不同尺度的低质量全色影像的第i行、第j列像素点的像素值；
Figure FDA0003412361170000031
为理想的高质量多光谱影像，即融合图像，
Figure FDA0003412361170000032
表示融合图像的g波段图像第i行、第j列像素点的像素值；
Figure FDA0003412361170000033
为低质量多光谱影像，其通过对原始多光谱影像进行采样得到，
Figure FDA0003412361170000034
表示低质量多光谱影像的g波段图像第i行、第j列像素点的像素值；i，j，g分别是行、列、波段的序号；全色影像与多光谱影像之间的数学关系描述为：
Pan＝k(g)MS(g)+b(g)T+E(g)，
其中，b(g)是多光谱影像第g波段对应的偏置系数，k(g)为多光谱影像第g波段对应的比例系数，T为维度与Pan相同、且其元素都为1的矩阵，E(g)是多光谱影像第g波段对应的残差矩阵，Pan是全色影像的像素值，MS(g)是g波段多光谱影像的像素值；通过低分辨率的全色多光谱影像拟合得到比例系数、偏置系数及残差矩阵，并利用这些系数得到融合的高分辨率全色影像。
2.如权利要求1所述的基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，其特征在于，所述的步骤S102，将步骤S101中计算得到的全色影像与多光谱影像的重叠区域中的全色影像读入CPU内存当中，以该重叠区域中的全色影像为基准，通过GPU对多光谱影像进行虚拟化，其具体步骤包括：
S1021，以全色影像的像方坐标为基准，将虚拟影像划分为均匀的规则格网，并且计算每个格网的四个角点的像方坐标；
S1022，对虚拟影像的每个格网的四个角点的像方坐标(s，l)，通过RFM模型计算其对应的物方坐标；
S1023，对虚拟影像的每个格网的四个角点的物方坐标，通过RFM模型计算多光谱影像中对应点的像方坐标(s′，l′)；
S1024，利用虚拟影像每个格网的四个角点的像方坐标，构建虚拟影像像方坐标与多光谱影像像方坐标的变换模型g()，其变换过程表示为(s，l)＝g(s′，l′)，该变换模型g()采用线性变换模型、仿射变换模型或透视变换模型中的一种；
S1025，对虚拟影像像素进行重采样；将虚拟影像像素重采样的过程映射到GPU进行处理，将GPU计算过程中的每一个CUDA线程设置为仅对一个像素进行计算；将虚拟影像上的每个像素点的坐标通过变换模型g()，计算得到多光谱影像上的像素点坐标，并且选用最近邻近采样、双线性采样或双三次重采样方法中的一种。
3.如权利要求1所述的基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，其特征在于，
所述的步骤S103，首先将虚拟化后的多光谱影像中的每个像素转换成与全色影像光学特性一致的虚拟全色单一波段像素，该过程表示为：
Mssp＝e1*BandB+e2*BandG+e3*BandR+e4*BandNIR，
式中，Mssp为转换后的与全色影像光学特性一致的虚拟全色单一波段像素值，e1、e2、e3、e4分别为蓝光波段B、绿光波段G、红光波段R、近红外波段NIR对应的波段范围占全色波段范围的比例，BandB、BandG、BandR、BandNIR分别代表多光谱虚拟影像中的蓝光波段、绿光波段、红光波段以及近红外波段影像的像素值；
提取全色影像与虚拟全色影像的同名点，该同名点作为多光谱虚拟影像进行数字微分纠正的输入；所述的提取全色影像与虚拟全色影像的同名点，采用相关系数法，以全色影像与虚拟全色影像的相关系数作为影像匹配的匹配测度，在指定的影像搜索范围中计算相关系数，而后选取相关系数最大的点作为同名点，相关系数是标准化的协方差函数；
选择M波段的全色影像上的某一像素点为初值点，以该初值点为中心，建立m×n个像素的匹配窗口，并在N波段的虚拟全色影像上建立与M波段的全色影像的初值点上同样大小和形状的匹配窗口，计算两个匹配窗口的相关系数，其计算公式为：
Figure FDA0003412361170000051
其中，R(c，r)表示两个匹配窗口在行和列上的坐标差分别为c和r时的相关系数，m，n分别表示匹配窗口的行和列大小，i，j分别表示匹配窗口内的行序号和列序号，gi，j表示M波段的全色影像和N波段的虚拟全色影像的匹配窗口内第i行、第j列的灰度值，gi+i，j+c表示上述的两个影像的匹配窗口内第i+r行、第j+c列的灰度值，c，r分别是N波段的虚拟全色影像的匹配窗口对于M波段的全色影像的匹配窗口中的同名点在行、列上的坐标差；对于已经进行了格网点匹配处理的多光谱影像，取c＝r＝0；对于未进行格网点匹配的影像，c和r的取值则根据相机的硬件参数值来确定；
在全色影像上从设定的格网点中确定控制点，采用相关系数法，依次获取每个格网点在虚拟全色影像上的同名点；将相关系数的计算过程映射到GPU上进行，GPU上的每个线程块对应一个格网点，线程块中的每个线程实现一组相关系数的计算，每个线程块计算完毕之后获取相关系数最大的同名点，并且利用抛物线法对相关系数最大的同名点相邻的四个点进行拟合，得到该相关系数最大的同名点的亚像素值。
4.如权利要求1所述的基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，其特征在于，
所述的步骤S104，对多光谱虚拟影像进行数字微分纠正，多光谱虚拟影像中的每个格网包括四组同名点，数字微分纠正过程f(x′，y′)的数学表达式为：
Figure FDA0003412361170000061
式中，(x，y)为所述的同名点在数字微分纠正前的像方坐标，(x′，y′)为所述的同名点在数字微分纠正后的像方坐标，g1、f1、g2、f2、g3、f3、g4、f4分别为数字微分纠正过程的8个变换系数，该8个变换系数通过对4组同名点联合建立的8个数字微分纠正过程方程进行解算得到；分别解算每个格网的变换系数，从而对多光谱虚拟影像进行数字微分纠正，再进行最近邻近采样、双线性采样或双三次重采样；将数字微分纠正的过程映射到GPU进行处理，将GPU计算过程中的每一个CUDA线程设置为仅对一个像素进行计算。
5.如权利要求1所述的基于CPU和GPU协同处理的全色和多光谱影像实时融合方法，其特征在于，步骤S105利用基于多尺度的平滑滤波亮度变换模型SFIM的融合算法实现；
所述的基于多尺度的SFIM模型的融合算法，是基于多尺度原理与SFIM模型的影像融合算法；SFIM模型的融合过程表示为，
Figure FDA0003412361170000062
其中，MS是多光谱影像的像素值，Pan是全色影像的像素值，Pan′是低分辨率全色影像像素值，Fusion是融合结果；Pan和Pan′的比值只保留了高分辨率图像的边缘细节信息；
对全色影像利用高斯金字塔方法，生成不同尺度的低分辨率全色影像，并向SFIM模型得到的融合影像中添加不同清晰程度的高频信息，以增强融合后的图像的清晰度。
CN202110877331.3A 2021-07-31 2021-07-31 基于cpu和gpu协同处理的全色和多光谱影像实时融合方法 Active CN113570536B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110877331.3A CN113570536B (zh) 2021-07-31 2021-07-31 基于cpu和gpu协同处理的全色和多光谱影像实时融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110877331.3A CN113570536B (zh) 2021-07-31 2021-07-31 基于cpu和gpu协同处理的全色和多光谱影像实时融合方法

Publications (2)

Publication Number Publication Date
CN113570536A CN113570536A (zh) 2021-10-29
CN113570536B true CN113570536B (zh) 2022-02-01

Family

ID=78169722

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110877331.3A Active CN113570536B (zh) 2021-07-31 2021-07-31 基于cpu和gpu协同处理的全色和多光谱影像实时融合方法

Country Status (1)

Country Link
CN (1) CN113570536B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114562982B (zh) * 2022-03-09 2023-09-26 北京市遥感信息研究所 一种光学和sar异源卫星影像联合平差的定权方法和装置
CN114418920B (zh) * 2022-03-30 2022-06-28 青岛大学附属医院 一种内窥镜多焦点图像融合方法
CN117078563B (zh) * 2023-10-16 2024-02-02 武汉大学 启明星一号卫星高光谱图像全色锐化方法及系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103049898A (zh) * 2013-01-27 2013-04-17 西安电子科技大学 带有薄云的多光谱和全色图像融合方法
CN103186893A (zh) * 2012-12-19 2013-07-03 中国科学院对地观测与数字地球科学中心 一种普适的高分辨率遥感图像融合方法
CN105160647A (zh) * 2015-10-28 2015-12-16 中国地质大学（武汉） 一种全色多光谱影像融合方法
CN106611410A (zh) * 2016-11-29 2017-05-03 北京空间机电研究所 基于金字塔模型的pansharpen融合优化方法
CN109118462A (zh) * 2018-07-16 2019-01-01 中国科学院东北地理与农业生态研究所 一种遥感影像融合方法
CN110930439A (zh) * 2019-12-04 2020-03-27 长光卫星技术有限公司 一种适用于高分辨率遥感影像高级产品自动生产系统

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108304775B (zh) * 2017-12-26 2021-06-11 北京市商汤科技开发有限公司 遥感图像识别方法、装置、存储介质以及电子设备
CN109544495B (zh) * 2018-11-13 2023-05-23 北京遥感设备研究所 一种基于高斯滤波和比值变换的SoC芯片图像处理融合方法
US10891527B2 (en) * 2019-03-19 2021-01-12 Mitsubishi Electric Research Laboratories, Inc. Systems and methods for multi-spectral image fusion using unrolled projected gradient descent and convolutinoal neural network
CN111681171B (zh) * 2020-06-15 2024-02-27 中国人民解放军军事科学院国防工程研究院 基于分块匹配的全色与多光谱图像高保真融合方法及装置

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103186893A (zh) * 2012-12-19 2013-07-03 中国科学院对地观测与数字地球科学中心 一种普适的高分辨率遥感图像融合方法
CN103049898A (zh) * 2013-01-27 2013-04-17 西安电子科技大学 带有薄云的多光谱和全色图像融合方法
CN105160647A (zh) * 2015-10-28 2015-12-16 中国地质大学（武汉） 一种全色多光谱影像融合方法
CN106611410A (zh) * 2016-11-29 2017-05-03 北京空间机电研究所 基于金字塔模型的pansharpen融合优化方法
CN109118462A (zh) * 2018-07-16 2019-01-01 中国科学院东北地理与农业生态研究所 一种遥感影像融合方法
CN110930439A (zh) * 2019-12-04 2020-03-27 长光卫星技术有限公司 一种适用于高分辨率遥感影像高级产品自动生产系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
CPU/GPU协同的光学卫星遥感数据高性能处理方法研究;方留杨;《中国优秀博硕士学位论文全文数据库(博士) 基础科学辑》;20160715;正文第2、3章 *
GPU Acceletated Automated feature extraction from satellite images;K Phani Tejaswi,D Shanmukha Rao;《International Journal of Distributed and Parallel Systems》;20130331;全文 *
基于稀疏表示的遥感图像融合;白晨帅，魏洁柔，苏维亚，王亚丁;《长春工程学院学报（自然科学版）》;20191230;全文 *

Also Published As

Publication number Publication date
CN113570536A (zh) 2021-10-29

Similar Documents

Publication Publication Date Title
CN113570536B (zh) 基于cpu和gpu协同处理的全色和多光谱影像实时融合方法
CN110415199B (zh) 基于残差学习的多光谱遥感图像融合方法及装置
CN109146787B (zh) 一种基于插值的双相机光谱成像系统的实时重建方法
CN110910437B (zh) 一种复杂室内场景的深度预测方法
CN104574347A (zh) 基于多源遥感数据的在轨卫星图像几何定位精度评价方法
CN106097252B (zh) 基于图Graph模型的高光谱图像超像素分割方法
CN111161199A (zh) 一种空谱融合的高光谱影像混合像元低秩稀疏分解方法
Liu et al. UMAG-Net: A new unsupervised multiattention-guided network for hyperspectral and multispectral image fusion
Zhang et al. Implicit neural representation learning for hyperspectral image super-resolution
CN114549669B (zh) 一种基于图像融合技术的彩色三维点云获取方法
CN111861880A (zh) 基于区域信息增强与块自注意力的图像超分与融合方法
CN116935214B (zh) 一种卫星多源遥感数据的时空谱融合方法
CN109741358B (zh) 基于自适应超图学习的超像素分割方法
CN113744136A (zh) 基于通道约束多特征融合的图像超分辨率重建方法和系统
CN112163990A (zh) 360度图像的显著性预测方法及系统
CN115760814A (zh) 一种基于双耦合深度神经网络的遥感图像融合方法及系统
CN111563866B (zh) 一种多源遥感图像融合方法
Liu et al. Circle-Net: An unsupervised lightweight-attention cyclic network for hyperspectral and multispectral image fusion
CN117474765B (zh) 基于参考影像纹理转移的dem超分辨率重建系统
CN112734636A (zh) 一种多源异构遥感影像的融合方法
CN109584194B (zh) 基于卷积变分概率模型的高光谱图像融合方法
Shen et al. An Overview of Image Super-resolution Reconstruction
Jin et al. Boosting single image super-resolution learnt from implicit multi-image prior
Qu et al. A principle design of registration-fusion consistency: Toward interpretable deep unregistered hyperspectral image fusion
Feng et al. The design of marine remote sensing image fusion model based on dynamic load balancing strategy

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