CN104794681B - 基于多冗余字典和稀疏重构的遥感图像融合方法 - Google Patents
基于多冗余字典和稀疏重构的遥感图像融合方法 Download PDFInfo
- Publication number
- CN104794681B CN104794681B CN201510208937.2A CN201510208937A CN104794681B CN 104794681 B CN104794681 B CN 104794681B CN 201510208937 A CN201510208937 A CN 201510208937A CN 104794681 B CN104794681 B CN 104794681B
- Authority
- CN
- China
- Prior art keywords
- mrow
- resolution
- mtd
- msub
- image block
- 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
- 238000007500 overflow downdraw method Methods 0.000 title claims abstract description 28
- 239000013598 vector Substances 0.000 claims description 58
- 239000011159 matrix material Substances 0.000 claims description 32
- 230000004927 fusion Effects 0.000 claims description 20
- 238000005070 sampling Methods 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 3
- 238000001228 spectrum Methods 0.000 abstract description 11
- 230000003595 spectral effect Effects 0.000 abstract description 7
- 238000012544 monitoring process Methods 0.000 abstract description 5
- 238000007499 fusion processing Methods 0.000 abstract description 3
- 230000007613 environmental effect Effects 0.000 abstract description 2
- 238000004611 spectroscopical analysis Methods 0.000 abstract 2
- 238000000034 method Methods 0.000 description 32
- 238000006243 chemical reaction Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- HUTDUHSNJYTCAR-UHFFFAOYSA-N ancymidol Chemical compound C1=CC(OC)=CC=C1C(O)(C=1C=NC=NC=1)C1CC1 HUTDUHSNJYTCAR-UHFFFAOYSA-N 0.000 description 2
- 239000004020 conductor Substances 0.000 description 2
- 230000007812 deficiency Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 241001269238 Data Species 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000003708 edge detection Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000002844 melting Methods 0.000 description 1
- 230000008018 melting Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4053—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4053—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution
- G06T3/4061—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution by injecting details from different spectral ranges
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于多冗余字典和稀疏重构的遥感图像融合方法，主要解决多光谱和全色图像融合过程中光谱信息和空间分辨率难以平衡的问题。其实现步骤为：1.对高分辨率全色图像和低分辨率多光谱图像分块；2.对低分辨率多光谱图像块的各个波段线性求和，得到低分辨率光谱‑强度图像块；3.根据低分辨率光谱‑强度图像块得到高分辨率光谱‑强度图像块；4.根据高分辨率光谱‑强度图像块和低分辨率多光谱图像块得到融合图像块；5.拼接融合图像块，得到融合图像。本发明具有融合图像的光谱保持性和空间分辨率高的优点，可用于到目标识别、气象监测、环境监测、城市规划以及防灾减灾。
Description
技术领域
本发明属于智能图像处理技术领域，特别涉及一种图像融合方法，可用于到目标识别、气象监测、环境监测、土地利用、城市规划以及防灾减灾。
背景技术
近年来，卫星技术在全球范围内迅猛发展，至今世界各主要发达国家和少数发展中国家，包括中国、印度等先后发射了数以百计的卫星，作业波段覆盖可见光至不可见的近红外、短波红外、中红外、远红外、微波等广阔频域。这些卫星每天向散布在世界各地的卫星地面站和移动接收站传送覆盖全球的海量卫星遥感数据。这些遥感数据为目标的识别、环境监测等提供了丰富而又宝贵的资料。但是由于实际应用中所需的遥感图像数据在时间、空间和光谱方面差异很大，而各种传感器提供的遥感图像数据都具有不同的特点，例如全色图像具有很高的空间分辨率，但是却没有丰富的光谱信息，而多光谱图像具有较好的光谱信息，但是分辨率却只有全色图像的四分之一，所以遥感技术应用的主要障碍不是数据源的不足，而是从这些数据源中提取更丰富、更有用和更可靠信息的能力大小。多源遥感图像融合是对图像数据之间的冗余性进行充分利用，从而降低在单一的遥感数据中存在的误差和不确定性，提高识别率和精确度。所以多源遥感数据的融合，尤其是多光谱和全色图像的融合，被认为是现代多源图像处理和分析中非常重要的一步。
目前，市场上使用的多光谱和全色图像融合方法主要有三类，一类是基于空间域的融合方法，一类是基于多尺度变换的融合方法，还有一类是目前很多学者在研究的基于压缩感知和字典学习的融合方法。
基于空间域的融合方法主要有HIS变换、PCA变换，Gram-Schmidt变换，以及Brovery变换等。这几种方法因其较低的计算复杂度，所以常被用于多个商业软件中。该类方法得到的融合图像具有较高的空间分辨率，但是光谱失真问题也是不可忽视的，而这个问题已经被很多的研究者所指出。
基于变换域的方法主要有基于拉普拉斯变换的融合方法、基于Wavelet变换的融合方法，以及基于多尺度几何分析的融合方法，如Contourlet、Bandlet和Shearlet等。这类方法相比于基于空间域的融合方法，融合图像的光谱失真情况得到改善，但是空间分辨率却没有基于空间域的融合方法好。这是由于融合图像的空间分辨率受到所采用的多尺度几何变换的约束，例如Wavelet变换只能将图像分解为三个方向，Shearlet变换等虽然能够对图像进行更多方向的分解，但还数目仍然是有限的。面对具有复杂细节信息的遥感图像，有限的方向很难达到最优的逼近，从而影响了融合图像的细节信息，造成融合图像空间分辨率的下降。
第三类是基于压缩感知和字典学习的融合方法。该类方法是在压缩感知的框架下，将多光谱和全色图像的融合问题转化为一个压缩感知重构问题。该类方法把低分辨率多光谱图像看做是融合图像的下采样图像，而把高分辨率全色图像看作是融合图像的各个波段的线性加和图像，从而将低分辨率多光谱图像和高分辨率全色图像作为压缩感知方法中的测量，然后设计相应的观测矩阵，通过求解一个优化问题，得到融合图像。该类方法的问题在于，当高分辨率全色图像近似等于融合图像的各个波段的线性加和的图像时，通过该方法能够得到较好的融合图像，而当高分辨率图像与融合图像的各个波段的线性加和的图像误差较大的时候，融合图像就会有较大的光谱失真。而在已有的文献中，该线性加和的权值都取为固定的数值，难以满足多样的遥感数据，造成融合图像光谱信息的丢失。因此对于多样的多光谱和全色图像，寻找一种更加有效的融合方法，是目前市场上急需解决的问题。
发明内容
本发明的目的在于针对上述已有技术的不足，提出一种基于多冗余字典和稀疏重构的遥感图像融合方法，以解决现有技术在多光谱和全色图像融合过程中光谱信息和空间分辨率难以平衡的问题，提高融合图像的质量。
为实现上述目的，本发明的实现步骤如下：
1)输入一幅高分辨率全色图像f和一幅低分辨率多光谱图像g，该低分辨率多光谱图像的大小是高分辨率全色图像大小的四分之一；
2)提取高分辨率全色图像f的Primal Sketch图，并根据该Primal Sketch图和几何模板，将高分辨率全色图像划为结构区域和非结构区域，并根据图像的方差统计特性，再将非结构区域划分为纹理区域和光滑区域，得到高分辨率全色图像的区域映射图；
3)对低分辨率多光谱图像g和高分辨率全色图像f进行分块，得到在相同位置的低分辨率多光谱图像的图像块gi,j和高分辨率全色图像的图像块记fj，其中i＝1,...,4,j＝1,...,J，i表示波段个数，j表示块的个数；
4)对低分辨率多光谱图像的图像块gi,j的各个波段进行线性求和，得到低分辨率光谱-强度图像块Ij，其中ωi表示低分辨率多光谱图像的第i个波段的权重值，m表示低分辨率多光谱图像的波段总个数；
5)根据低分辨率光谱-强度图像块Ij和高分辨率全色图像块fj，得到高分辨率光谱-强度图像块列向量J'j；
6)根据高分辨率光谱-强度图像的列向量J'j和低分辨率多光谱图像块gi,j，得到融合图像块Pi,j：
6a)将步骤2)得到的区域映射图映射到低分辨率多光谱图像上，并进行判断：如果低分辨率多光谱图像块gi,j在光滑区域，则执行步骤6b)，否则，执行步骤6c)；
6b)对低分辨率多光谱图像块gi,j进行上采样得到融合图像块Pi,j，其中i表示融合图像块的第i个波段，j表示融合图像的第j块；
6c)根据低分辨率多光谱图像块gi,j组成低分辨率多光谱图像块列向量其中g'i,j是低分辨率多光谱图像块gi,j的列向量；
6d)基于高分辨率光谱-强度图像块列向量J'j和低分辨率多光谱图像块列向量g'j，通过BP算法求解式，得到融合图像块列向量Pj'的稀疏系数
其中，表示·的二范数的平方,Z为融合图像块重构字典，其根据低分辨率多光谱图像块gi,j所在区域选择，λ2是稀疏项的参数,β2是约束项的参数，δj表示融合图像块在融合图像块重构字典Z下的稀疏系数，M2是融合图像块列向量的下采样矩阵，M3是权值矩阵：M3＝(ω1Kω2Lω3Lω4L)，
式中F2×2是2×2的单位矩阵，E8×8是8×8的单位矩阵，L是4行1列的矩阵，矩阵中的值都是1，表示内积运算，K是64×64的单位矩阵，ω1,ω2,ω3,ω4分别是第一波段、第二波段、第三波段和第四波段的权值，ω1,ω2,ω3,ω4的值均取为0.25；
6e)根据稀疏系数和融合图像块重构字典Z，计算得到融合图像块的列向量Pj'，并将其转成融合图像的图像块Pi,j；
7)将所有的融合图像块Pi,j拼接，得到融合图像。
本发明与现有技术相比，具有如下效果：
(a)本发明由于使用多冗余字典和稀疏重构的方法对低分辨率多光谱和高分辨率全色图像进行融合，将图像融合的问题转化为一个稀疏重构问题，克服了传统的融合方法得到的融合图像中空间分辨率和光谱信息难以平衡的问题，提高了融合图像的质量。
(b)本发明对融合图像块重构字典的选择，是根据低分辨率多光谱图像块所在区域选择，对不同的区域选择不同的固定基冗余字典对融合图像进行重构，克服了传统的融合图像重构方法中采用单一的固定基字典不能较好的描述各个区域的特征，以及传统的字典学习方法中训练样本难以直接获取的问题，从而提高了融合图像的空间分辨率。
附图说明
图1是本发明的多光谱图像和全色图像融合流程框图；
图2是本发明中将高分辨率全色图像划分为结构与非结构区的示意图；
图3是本发明中对不同区域设计的部分字典结构图；
图4是用本发明和现有方法中对农田图像的融合结果；
图5是用本发明和现有技术中对城区图像的融合结果。
具体实施方式
以下参照附图，对本发明的技术方案和效果做进一步详细描述。
参照图1，本发明的实现步骤如下：
步骤1，输入图像。
从QuickBird卫星图像中选择一幅高分辨率全色图像f和一幅大小是高分辨率全色图像大小的四分之一低分辨率多光谱图像g，作为输入图像。
步骤2，获得高分辨率全色图像的区域映射图。
2a)划分高分辨率全色图像的结构与非结构区域
区域的划分可以根据传统的边缘检测方法、区域生长，和Primal Sketch图等方法得到，本发明根据Primal Sketch图得到，即先从高分辨率全色图像先提取高分辨率全色图像f的Primal Sketch图，Primal Sketch图如图2(b)所示，该图中的每一条线段都是有方向信息的，包含了图像中的边、线结构，由于Primal Sketch图中对于边、线结构的描述都是由一个像素表示的，而图像的边缘具有连续性，其由3到5个像素描述，所以根据PrimalSketch图中线段的方向，以线段上的点为中心，沿着该条线段的方向，构造一个7×7的几何模版；再根据该几何模版，对Primal Sketch图中线段上的每一点都应用该模版，得到高分辨率全色图像的结构与非结构区域，如图2(c)所示；
2b)根据图像的方差统计特性，将非结构区域划分为纹理区域和光滑区域，得到高分辨率全色图像的区域映射图：
2b1)以高分辨率全色图像的非结构区域的每一个像素点(i',j')为中心，取5×5的窗口，计算窗口内像素点的方差值V(i',j')：
其中f(p,q)表示高分辨率全色图像f在(p,q)点的值，是以像素点(i',j')为中心的5×5窗口内的像素值得均值，
2b2)设定阈值T＝2，对高分辨率全色图像非结构区域的所有像素点进行判定:如果V(i',j')＞T，则像素点(i',j')被判定在纹理区域；否则，被判定在光滑区域，从而得到高分辨率全色图像的区域映射图。
步骤3，分别对低分辨率多光谱图像g和高分辨率全色图像f进行分块。
分块是通过对整个图像滑窗得到的：
对低分辨率多光谱图像g，其滑窗窗口的大小为2×2，相邻窗口的步长为1，得到的低分辨率多光谱图像的图像块为gi,j，i＝1,...,4,j＝1,...,J，i表示波段个数，j表示块的个数；
对高分辨率全色图像f，其滑窗窗口的大小为8×8，相邻窗口的步长为4，得到高分辨率全色图像的图像块fj。
步骤4，根据低分辨率多光谱图像的图像块gi,j得到低分辨率光谱-强度图像块Ij。
对低分辨率多光谱图像的图像块gi,j的各个波段进行线性求和，得到低分辨率光谱-强度图像块：其中ωi表示低分辨率多光谱图像的第i个波段的权重值，ωi＝0.25。
步骤5，根据低分辨率光谱-强度图像块Ij和高分辨率全色图像块fj，得到高分辨率光谱-强度图像块列向量J'j。
将步骤2得到的区域映射图映射到低分辨率光谱-强度图像上，并进行判断：
如果低分辨率光谱-强度图像块Ij在光滑区域，则对低分辨率光谱-强度图像块Ij进行上采样，得到高分辨率光谱-强度图像块Jj，并将其转成高分辨率光谱-强度图像块列向量J'j，其中采样的方法有最近邻插值，线性插值，双线性插值，立方插值，本实例采用线性插值法；
如果低分辨率光谱-强度图像块Ij在纹理或结构区域，则执行如下步骤：
5a)分别将低分辨率光谱-强度图像块Ij和高分辨率全色图像块fj转成低分辨率光谱-强度图像块列向量Ij'和高分辨率全色图像块列向量fj'，并基于这两个列向量，通过BP算法求解式，得到高分辨率光谱-强度图像块列向量J'j的稀疏系数
其中表示·的二范数的平方,αj表示高分辨率光谱-强度图像块列向量J'j在高分辨率光谱-强度图像块重构字典D下的稀疏系数,参数λ1＝1，β1＝0.01，D是高分辨率光谱-强度图像块重构字典，λ1稀疏项的参数,β1是约束项的参数，M1是高分辨率光谱-强度图像块列向量的下采样矩阵，MT是梯度矩阵，
式中F2×2是2×2的单位矩阵，L是4行1列的矩阵，矩阵中的值都是1，表示内积运算，M表示高分辨率光谱-强度图像块列向量J'j的高；
所述高分辨率光谱-强度图像块重构字典D，是根据低分辨率光谱-强度图像块Ij所在区域选择：
如果低分辨率光谱-强度图像块Ij在纹理区域，则高分辨率光谱-强度图像块重构字典D选择为DCT字典D1，部分DCT字典D1如图3(a)所示，并根据该字典计算得到稀疏系数
如果低分辨率光谱-强度图像块Ij在结构区域，则高分辨率光谱-强度图像块重构字典D选择Curvelet字典D2和Ridgelet字典D3，并分别根据Curvelet字典D2和Ridgelet字典D3计算得到稀疏系数和部分Curvelet字典D2如图3(b)所示，部分Ridgelet字典D3如图3(c)所示；
5b)根据稀疏系数和高分辨率光谱-强度图像块重构字典D，计算得到高分辨率光谱-强度图像块的列向量J'j：
对选择的高分辨率光谱-强度图像块重构字典D进行判断：
如果高分辨率光谱-强度图像块重构字典D选择DCT字典D1，则计算高分辨率光谱-强度图像块列向量
如果高分辨率光谱-强度图像块重构字典D选择Curvelet字典D2和Ridgelet字典D3，则根据Curvelet字典D2和Ridgelet字典D3分别计算得到两个列向量和并根据下式得到高分辨率光谱-强度图像块的列向量J'j：
步骤6，根据高分辨率光谱-强度图像的图像块列向量J'j和低分辨率多光谱图像块gi,j，得到融合图像块Pi,j。
将步骤2得到的区域映射图映射到低分辨率多光谱图像上，并进行判断：
如果低分辨率多光谱图像块gi,j在光滑区域，则对低分辨率多光谱图像块gi,j进行上采样，到融合图像块Pi,j，所述上采样的方法有最近邻插值，线性插值，双线性插值，立方插值，本实例采用线性插值；
如果低分辨率多光谱图像块gi,j在纹理和结构区域，则执行如下步骤：
6a)根据低分辨率多光谱图像块gi,j组成低分辨率多光谱图像块列向量：其中g'i,j是低分辨率多光谱图像块gi,j的列向量；
6b)基于高分辨率光谱-强度图像块列向量J'j和低分辨率多光谱图像块列向量g'j，通过BP算法求解式，得到融合图像块列向量Pj'的稀疏系数
其中，表示·的二范数的平方,Z是融合图像块重构字典，λ2是稀疏项的参数,β2是约束项的参数，M2是融合图像块列向量的下采样矩阵，M3是权值矩阵：M3＝(ω1Kω2Lω3Lω4L)，
式中F2×2是2×2的单位矩阵，E8×8是8×8的单位矩阵，L是4行1列的矩阵，矩阵中的值都是1，表示内积运算，K是64×64的单位矩阵；
所述融合图像块重构字典Z，是根据低分辨率多光谱图像块gi,j所在区域选择：
如果低分辨率多光谱图像块gi,j在纹理区域，则融合图像块重构字典Z选择为多冗余DCT字典并根据该字典计算得到稀疏系数
如果低分辨率多光谱图像块gi，j在结构区域，则融合图像块重构字典Z选择为多冗余Curvelet字典和多冗余Ridgelet字典并分别根据多冗余Curvelet字典Z2和多冗余Ridgelet字典Z3计算得到稀疏系数和
6c)根据稀疏系数和融合图像块重构字典Z，计算得到融合图像块的列向量Pj'：
对融合图像块重构字典Z进行判断：
如果融合图像块重构字典Z选择选择为多冗余DCT字典Z1，则计算融合图像块列向量
如果融合图像块重构字典Z选择为多冗余Curvelet字典Z2和多冗余Ridgelet字典Z3，则根据多冗余Curvelet字典Z2和多冗余Ridgelet字典Z3分别计算得到两个列向量和并根据下式得到融合图像块的列向量Pj'：
6d)将融合图像块的列向量Pj'转成融合图像的图像块Pi,j。
步骤7，将所有的融合图像块Pi,j进行拼接，得到最终的融合图像。
本发明的效果可以用下列的仿真实验进一步说明:
(1)仿真条件
本发明的仿真的硬件条件为：windows XP,SPI,CPU Pentium(R)4,基本频率2.4GHZ,软件平台为：MatlabR2012a，仿真选用的图片来源是QuickBrid卫星图像的低分辨率多光谱和高分辨率全色图像，分别对应图4和图5，其中图4(a)和图5(a)是QuickBrid卫星图像的低分辨率多光谱图像，图4(b)和图5(b)是QuickBrid卫星图像的高分辨率全色图像。
仿真方法分别用本发明方法和GS的融合方法、Brovery融合方法、P+XS融合方法，DWT融合方法对低分辨率多光谱图像和高分辨率全色图像进行融合。
(2)仿真内容及结果
仿真1，用本发明和现有的四种方法对图4(a)和图4(b)进行仿真，得到如图4(c)-(g)所示的融合图像。其中图4(c)是GS法的融合结果图，图4(d)是Brovery法的融合结果图，图4(e)是P+XS法的融合结果图，图4(f)是DWT法的融合结果图，图4(g)是本发明融合的结果图，从图4(c)-(g)的融合结果图可见，本发明融合的结果图更清晰，图像的光谱信息更好。
仿真2，用本发明和现有的四种方法对图5(a)和图5(b)进行仿真，得到如图5(c)-(g)所示的融合图像。其中图5(c)是GS法的融合结果图，图5(d)是Brovery法的融合结果图，图5(e)是P+XS法的融合结果图，图5(f)是DWT法的融合结果图，图5(g)是本发明融合的结果图，从图5(c)-(g)的融合结果图可见，本发明融合的结果图更清晰，图像的光谱信息更好。
以上实验结果表明：本发明相比现有技术在解决多光谱和全色图像融合过程中光谱信息和空间分辨率难以平衡，及多光谱和全色图像融合后光谱失真或者空间分辨率不高的问题上，具有融合图像的光谱保持性和空间分辨率高的优点，提高了融合图像的质量。
Claims (8)
1.一种基于多冗余字典和稀疏重构的遥感图像融合方法，包括如下步骤：
1)输入一幅高分辨率全色图像f和一幅低分辨率多光谱图像g，该低分辨率多光谱图像的大小是高分辨率全色图像大小的四分之一；
2)提取高分辨率全色图像f的Primal Sketch图，并根据该Primal Sketch图和几何模板，将高分辨率全色图像划为结构区域和非结构区域，并根据图像的方差统计特性，再将非结构区域划分为纹理区域和光滑区域，得到高分辨率全色图像的区域映射图；
3)对低分辨率多光谱图像g和高分辨率全色图像f进行分块，得到在相同位置的低分辨率多光谱图像的图像块gi,j和高分辨率全色图像的图像块记fj，其中i＝1,...,4,j＝1,...,J，i表示波段个数，j表示块的个数；
4)对低分辨率多光谱图像的图像块gi,j的各个波段进行线性求和，得到低分辨率光谱-强度图像块Ij，其中ωi表示低分辨率多光谱图像的第i个波段的权重值，m表示低分辨率多光谱图像的波段总个数：
5)根据低分辨率光谱-强度图像块Ij和高分辨率全色图像块fj，得到高分辨率光谱-强度图像块列向量J'j；
6)根据高分辨率光谱-强度图像的列向量J'j和低分辨率多光谱图像块gi,j，得到融合图像块Pi,j：
6a)将步骤2)得到的区域映射图映射到低分辨率多光谱图像上，并进行判断：如果低分辨率多光谱图像块gi,j在光滑区域，则执行步骤6b)，否则，执行步骤6c)；
6b)对低分辨率多光谱图像块gi,j进行上采样得到融合图像块Pi,j，其中i表示融合图像块的第i个波段，j表示融合图像的第j块；
6c)根据低分辨率多光谱图像块gi,j组成低分辨率多光谱图像块列向量其中g'i,j是低分辨率多光谱图像块gi,j的列向量；
6d)基于高分辨率光谱-强度图像块列向量J'j和低分辨率多光谱图像块列向量g'j，通过BP算法求解式，得到融合图像块列向量Pj'的稀疏系数
<mrow>
<msub>
<mover>
<mi>&delta;</mi>
<mo>^</mo>
</mover>
<mi>j</mi>
</msub>
<mo>=</mo>
<munder>
<mrow>
<mi>arg</mi>
<mi>min</mi>
</mrow>
<msub>
<mi>&delta;</mi>
<mi>j</mi>
</msub>
</munder>
<msub>
<mi>&lambda;</mi>
<mn>2</mn>
</msub>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>&delta;</mi>
<mi>j</mi>
</msub>
<mo>|</mo>
<msub>
<mo>|</mo>
<mn>1</mn>
</msub>
<mo>+</mo>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mo>|</mo>
<msubsup>
<mi>g</mi>
<mi>j</mi>
<mo>&prime;</mo>
</msubsup>
<mo>-</mo>
<msub>
<mi>M</mi>
<mn>2</mn>
</msub>
<msub>
<mi>Z&delta;</mi>
<mi>j</mi>
</msub>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msub>
<mi>&beta;</mi>
<mn>2</mn>
</msub>
<mo>|</mo>
<mo>|</mo>
<msubsup>
<mi>J</mi>
<mi>j</mi>
<mo>&prime;</mo>
</msubsup>
<mo>-</mo>
<msub>
<mi>M</mi>
<mn>3</mn>
</msub>
<msub>
<mi>Z&delta;</mi>
<mi>j</mi>
</msub>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
其中，表示·的二范数的平方,Z为融合图像块重构字典，其根据低分辨率多光谱图像块gi,j所在区域选择，λ2是稀疏项的参数,β2是约束项的参数，M2是融合图像块列向量的下采样矩阵，M3是权值矩阵：M3＝(ω1Kω2Lω3Lω4L)，δj表示融合图像块在融合图像块重构字典Z下的稀疏系数；
式中F2×2是2×2的单位矩阵，E8×8是8×8的单位矩阵，L是4行1列的矩阵，矩阵中的值都是1，表示内积运算，K是64×64的单位矩阵，ω1,ω2,ω3,ω4分别是第一波段、第二波段、第三波段和第四波段的权值，ω1,ω2,ω3,ω4的值均取为0.25；
6e)根据稀疏系数和融合图像块重构字典Z，计算得到融合图像块的列向量Pj'，并将其转成融合图像的图像块Pi,j；
7)将所有的融合图像块Pi,j拼接，得到融合图像。
2.根据权利要求1所述的基于多冗余字典和稀疏重构的遥感图像融合方法，其中所述步骤2)中将非结构区域划分为纹理区域和光滑区域，按照如下步骤进行：
2a)以高分辨率全色图像的非结构区域的每一个像素点(i',j')为中心，取5×5的窗口，计算窗口内像素点的方差值V(i',j')：
<mrow>
<mi>V</mi>
<mrow>
<mo>(</mo>
<msup>
<mi>i</mi>
<mo>&prime;</mo>
</msup>
<mo>,</mo>
<msup>
<mi>j</mi>
<mo>&prime;</mo>
</msup>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>5</mn>
<mo>&times;</mo>
<mn>5</mn>
</mrow>
</mfrac>
<msqrt>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>p</mi>
<mo>=</mo>
<msup>
<mi>i</mi>
<mo>&prime;</mo>
</msup>
<mo>-</mo>
<mn>2</mn>
</mrow>
<mrow>
<msup>
<mi>i</mi>
<mo>&prime;</mo>
</msup>
<mo>+</mo>
<mn>2</mn>
</mrow>
</munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>p</mi>
<mo>=</mo>
<msup>
<mi>j</mi>
<mo>&prime;</mo>
</msup>
<mo>-</mo>
<mn>2</mn>
</mrow>
<mrow>
<msup>
<mi>j</mi>
<mo>&prime;</mo>
</msup>
<mo>+</mo>
<mn>2</mn>
</mrow>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mi>f</mi>
<mo>(</mo>
<mrow>
<mi>p</mi>
<mo>,</mo>
<mi>q</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<mover>
<mi>F</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mrow>
其中f(p,q)表示高分辨率全色图像f在(p,q)点的值，是以像素点(i',j')为中心的5×5窗口内的像素值得均值，
2b)设定阈值T＝2，对高分辨率全色图像非结构区域的所有像素点进行判定:
如果V(i',j')＞T，则像素点(i',j')被判定在纹理区域；否则，被判定在光滑区域。
3.根据权利要求1所述的基于多冗余字典和稀疏重构的遥感图像融合方法，
其中所述步骤3)中对低分辨率多光谱图像g和高分辨率全色图像f进行分块，是通过对整个图像滑窗得到的，即对低分辨率多光谱图像g，滑窗窗口的大小为2×2，相邻窗口的步长为1，得到的低分辨率多光谱图像的图像块gi,j；对高分辨率全色图像f，滑窗窗口的大小为8×8，相邻窗口的步长为4，得到高分辨率全色图像的图像块fj。
4.根据权利要求1所述的基于多冗余字典和稀疏重构的遥感图像融合方法，其中所述步骤5)，按如下步骤实现：
5a)将步骤1和步骤2)得到的区域映射图映射到低分辨率光谱-强度图像上，并进行判断：如果低分辨率光谱-强度图像块Ij在光滑区域，则执行步骤5b)，否则，执行步骤5c)；
5b)对低分辨率光谱-强度图像块Ij上采样得到高分辨率光谱-强度图像块Jj，并将其转成高分辨率光谱-强度图像块列向量J'j；
5c)分别将低分辨率光谱-强度图像块Ij和高分辨率全色图像块fj转成低分辨率光谱-强度图像块列向量Ij'和高分辨率全色图像块列向量fj'，并基于这两个列向量，通过BP算法求解式，得到高分辨率光谱-强度图像块列向量J'j的稀疏系数
<mrow>
<msub>
<mover>
<mi>&alpha;</mi>
<mo>^</mo>
</mover>
<mi>j</mi>
</msub>
<mo>=</mo>
<munder>
<mrow>
<mi>arg</mi>
<mi>min</mi>
</mrow>
<msub>
<mi>&alpha;</mi>
<mi>j</mi>
</msub>
</munder>
<msub>
<mi>&lambda;</mi>
<mn>1</mn>
</msub>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>&alpha;</mi>
<mi>j</mi>
</msub>
<mo>|</mo>
<msub>
<mo>|</mo>
<mn>1</mn>
</msub>
<mo>+</mo>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mo>|</mo>
<msup>
<msub>
<mi>I</mi>
<mi>j</mi>
</msub>
<mo>&prime;</mo>
</msup>
<mo>-</mo>
<msub>
<mi>M</mi>
<mn>1</mn>
</msub>
<msub>
<mi>D&alpha;</mi>
<mi>j</mi>
</msub>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msub>
<mi>&beta;</mi>
<mn>1</mn>
</msub>
<mo>|</mo>
<mo>|</mo>
<msub>
<mi>M</mi>
<mi>T</mi>
</msub>
<msubsup>
<mi>f</mi>
<mi>j</mi>
<mo>&prime;</mo>
</msubsup>
<mo>-</mo>
<msub>
<mi>M</mi>
<mi>T</mi>
</msub>
<msub>
<mi>D&alpha;</mi>
<mi>j</mi>
</msub>
<mo>|</mo>
<msubsup>
<mo>|</mo>
<mn>2</mn>
<mn>2</mn>
</msubsup>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
其中表示·的二范数的平方,αj表示高分辨率光谱-强度图像块列向量J'j在高分辨率光谱-强度图像块重构字典D下的稀疏系数,D是高分辨率光谱-强度图像块重构字典，λ1稀疏项的参数,β1是约束项的参数，M1是高分辨率光谱-强度图像块列向量的下采样矩阵，MT是梯度矩阵，
<mrow>
<msub>
<mi>M</mi>
<mi>T</mi>
</msub>
<mo>=</mo>
<msub>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mn>2</mn>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mn>2</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow></mrow>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mrow>
<mi>M</mi>
<mo>&times;</mo>
<mi>M</mi>
</mrow>
</msub>
<mo>,</mo>
</mrow>
式中F2×2是2×2的单位矩阵，L是4行1列的矩阵，矩阵中的值都是1，表示内积运算，M表示高分辨率光谱-强度图像块列向量J'j的高；
5d)根据稀疏系数和高分辨率光谱-强度图像块重构字典D，计算得到高分辨率光谱-强度图像块的列向量J'j。
5.根据权利要求4所述的基于多冗余字典和稀疏重构的遥感图像融合方法，其中所述步骤5c)中的高分辨率光谱-强度图像块重构字典D，是根据低分辨率光谱-强度图像块Ij所在区域选择：
如果低分辨率光谱-强度图像块Ij在纹理区域，则高分辨率光谱-强度图像块重构字典D选择为DCT字典D1，并根据该字典计算得到稀疏系数
如果低分辨率光谱-强度图像块Ij在结构区域，则高分辨率光谱-强度图像块重构字典D选择Curvelet字典D2和Ridgelet字典D3，并分别根据Curvelet字典D2和Ridgelet字典D3计算得到稀疏系数和
6.根据权利要求4或5所述的基于多冗余字典和稀疏重构的遥感图像融合方法，其中所述步骤5d)中根据稀疏系数和高分辨率光谱-强度图像块重构字典D，计算得到高分辨率光谱-强度图像块列向量J'j，按照如下步骤进行：
5d1)根据高分辨率光谱-强度图像块重构字典D进行判断：如果高分辨率光谱-强度图像块重构字典D选择DCT字典D1，则执行步骤5d2)；如果高分辨率光谱-强度图像块重构字典D选择Curvelet字典D2和Ridgelet字典D3，执行步骤5d3)；
5d2)计算高分辨率光谱-强度图像块列向量
5d3)根据Curvelet字典D2和Ridgelet字典D3分别计算得到两个列向量和并根据下式得到高分辨率光谱-强度图像块的列向量J'j：
7.根据权利要求1所述的基于多冗余字典和稀疏重构的遥感图像融合方法，其中所述步骤6d)中根据低分辨率多光谱图像块gi,j所在区域选择融合图像块重构字典Z，按如下步骤进行：
如果低分辨率多光谱图像块gi,j在纹理区域，则融合图像块重构字典Z选择为多冗余DCT字典并根据该字典计算得到稀疏系数
如果低分辨率多光谱图像块gi,j在结构区域，则融合图像块重构字典Z选择为多冗余Curvelet字典和多冗余Ridgelet字典并分别根据多冗余Curvelet字典Z2和多冗余Ridgelet字典Z3计算得到稀疏系数和
8.根据权利要求1或7所述的基于多冗余字典和稀疏重构的遥感图像融合方法，其中所述步骤6e)中根据稀疏系数和融合图像块重构字典Z得到融合图像块列向量Pj'，按照如下步骤进行：
6e1)根据融合图像块重构字典Z判断：如果融合图像块重构字典Z选择选择为多冗余DCT字典Z1，则执行步骤6e2)；如果融合图像块重构字典Z选择为多冗余Curvelet字典Z2和多冗余Ridgelet字典Z3，执行步骤6e3)；
6e2)计算融合图像块列向量
6e3)根据多冗余Curvelet字典Z2和多冗余Ridgelet字典Z3分别计算得到两个列向量和并根据下式得到融合图像块的列向量Pj'：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510208937.2A CN104794681B (zh) | 2015-04-28 | 2015-04-28 | 基于多冗余字典和稀疏重构的遥感图像融合方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510208937.2A CN104794681B (zh) | 2015-04-28 | 2015-04-28 | 基于多冗余字典和稀疏重构的遥感图像融合方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104794681A CN104794681A (zh) | 2015-07-22 |
| CN104794681B true CN104794681B (zh) | 2018-03-13 |
Family
ID=53559464
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510208937.2A Active CN104794681B (zh) | 2015-04-28 | 2015-04-28 | 基于多冗余字典和稀疏重构的遥感图像融合方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104794681B (zh) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105701769A (zh) * | 2016-01-13 | 2016-06-22 | 河海大学 | 边界灰度分布相关的合成孔径雷达遥感图像分块重构方法 |
| CN105719262B (zh) * | 2016-01-21 | 2018-06-22 | 西北大学 | 基于子字典稀疏重构的全色与多光谱遥感图像融合方法 |
| CN106251320B (zh) * | 2016-08-15 | 2019-03-26 | 西北大学 | 基于联合稀疏与结构字典的遥感图像融合方法 |
| CN107169945B (zh) * | 2017-04-25 | 2019-06-21 | 西安电子科技大学 | 基于稀疏张量和多视图特征的遥感图像融合方法 |
| CN107705280B (zh) * | 2017-10-23 | 2020-12-15 | 北京航空航天大学 | 一种结构驱动的光谱映射遥感图像融合方法 |
| CN109102480B (zh) * | 2018-07-06 | 2022-02-22 | 中科星图股份有限公司 | 一种适用于分布式架构的Gram-Schmidt融合方法 |
| CN111192196A (zh) * | 2019-12-23 | 2020-05-22 | 中电健康云科技有限公司 | 一种提高推扫式光谱仪高光谱图像的实时分辨率的方法 |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101996396A (zh) * | 2010-09-16 | 2011-03-30 | 湖南大学 | 一种基于压缩传感理论的卫星遥感图像融合方法 |
| CN102013106A (zh) * | 2010-10-27 | 2011-04-13 | 西安电子科技大学 | 基于Curvelet冗余字典的图像稀疏表示方法 |
| CN102542549A (zh) * | 2012-01-04 | 2012-07-04 | 西安电子科技大学 | 基于压缩感知的多光谱与全色图像超分辨融合方法 |
| CN102609910A (zh) * | 2012-01-13 | 2012-07-25 | 西安电子科技大学 | 基于Ridgelet冗余字典的遗传进化图像重构方法 |
| CN102651124A (zh) * | 2012-04-07 | 2012-08-29 | 西安电子科技大学 | 基于冗余字典稀疏表示和评价指标的图像融合方法 |
| CN103208102A (zh) * | 2013-03-29 | 2013-07-17 | 上海交通大学 | 一种基于稀疏表示的遥感图像融合方法 |
| CN103295198A (zh) * | 2013-05-13 | 2013-09-11 | 西安电子科技大学 | 基于冗余字典和结构稀疏的非凸压缩感知图像重构方法 |
| GB2501171A (en) * | 2012-03-12 | 2013-10-16 | Toshiba Kk | A parallel multi-frame superresolution image processing system |
| CN103914817A (zh) * | 2014-03-04 | 2014-07-09 | 西安电子科技大学 | 一种基于区域划分和插值的多光谱和全色图像融合方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040151376A1 (en) * | 2003-02-05 | 2004-08-05 | Konica Minolta Holdings, Inc. | Image processing method, image processing apparatus and image processing program |
-
2015
- 2015-04-28 CN CN201510208937.2A patent/CN104794681B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101996396A (zh) * | 2010-09-16 | 2011-03-30 | 湖南大学 | 一种基于压缩传感理论的卫星遥感图像融合方法 |
| CN102013106A (zh) * | 2010-10-27 | 2011-04-13 | 西安电子科技大学 | 基于Curvelet冗余字典的图像稀疏表示方法 |
| CN102542549A (zh) * | 2012-01-04 | 2012-07-04 | 西安电子科技大学 | 基于压缩感知的多光谱与全色图像超分辨融合方法 |
| CN102609910A (zh) * | 2012-01-13 | 2012-07-25 | 西安电子科技大学 | 基于Ridgelet冗余字典的遗传进化图像重构方法 |
| GB2501171A (en) * | 2012-03-12 | 2013-10-16 | Toshiba Kk | A parallel multi-frame superresolution image processing system |
| CN102651124A (zh) * | 2012-04-07 | 2012-08-29 | 西安电子科技大学 | 基于冗余字典稀疏表示和评价指标的图像融合方法 |
| CN103208102A (zh) * | 2013-03-29 | 2013-07-17 | 上海交通大学 | 一种基于稀疏表示的遥感图像融合方法 |
| CN103295198A (zh) * | 2013-05-13 | 2013-09-11 | 西安电子科技大学 | 基于冗余字典和结构稀疏的非凸压缩感知图像重构方法 |
| CN103914817A (zh) * | 2014-03-04 | 2014-07-09 | 西安电子科技大学 | 一种基于区域划分和插值的多光谱和全色图像融合方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104794681A (zh) | 2015-07-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104794681B (zh) | 基于多冗余字典和稀疏重构的遥感图像融合方法 | |
| US11521377B1 (en) | Landslide recognition method based on laplacian pyramid remote sensing image fusion | |
| Ranchin et al. | Image fusion—The ARSIS concept and some successful implementation schemes | |
| CN113128134B (zh) | 一种矿区生态环境演变驱动因子权重量化分析方法 | |
| CN103218796B (zh) | 一种全色—多光谱遥感图像融合方法 | |
| CN113420662A (zh) | 基于孪生多尺度差异特征融合的遥感影像变化检测方法 | |
| CN101996396A (zh) | 一种基于压缩传感理论的卫星遥感图像融合方法 | |
| CN111738111A (zh) | 基于多分支级联空洞空间金字塔的高分辨遥感图像的道路提取方法 | |
| CN107392130A (zh) | 基于阈值自适应和卷积神经网络的多光谱图像分类方法 | |
| CN112819737B (zh) | 基于3d卷积的多尺度注意力深度卷积网络的遥感图像融合方法 | |
| CN106384332A (zh) | 基于Gram‑Schmidt的无人机影像与多光谱影像融合方法 | |
| CN114565843A (zh) | 一种时间序列遥感图像融合方法 | |
| CN107563411A (zh) | 基于深度学习的在线sar目标检测方法 | |
| CN108932710A (zh) | 遥感时空信息融合方法 | |
| Liu et al. | A shallow-to-deep feature fusion network for VHR remote sensing image classification | |
| CN112017160A (zh) | 一种基于多策略组合的多源遥感影像道路材质精细提取方法 | |
| CN103440500A (zh) | 高光谱遥感图像分类与识别方法 | |
| CN116258976A (zh) | 一种分层次Transformer的高分辨率遥感图像语义分割方法及系统 | |
| CN113625363A (zh) | 伟晶岩型锂矿的找矿方法及装置、计算机设备及介质 | |
| CN115457396B (zh) | 一种基于遥感影像的地表目标地物检测方法 | |
| CN116309070A (zh) | 一种高光谱遥感图像超分辨率重建方法、装置及计算机设备 | |
| CN116091833A (zh) | 注意力与Transformer高光谱图像分类方法及系统 | |
| CN116029904A (zh) | 一种基于LRepSR与I-E知识蒸馏的遥感图像超分辨率重建方法 | |
| CN112529828B (zh) | 参考数据非敏感的遥感影像时空融合模型构建方法 | |
| CN103914817B (zh) | 一种基于区域划分和插值的多光谱和全色图像融合方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| EXSB | Decision made by sipo to initiate substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |