CN114926452B - 一种基于NSST和β散度非负矩阵分解的遥感图像融合方法 - Google Patents
一种基于NSST和β散度非负矩阵分解的遥感图像融合方法 Download PDFInfo
- Publication number
- CN114926452B CN114926452B CN202210657229.7A CN202210657229A CN114926452B CN 114926452 B CN114926452 B CN 114926452B CN 202210657229 A CN202210657229 A CN 202210657229A CN 114926452 B CN114926452 B CN 114926452B
- Authority
- CN
- China
- Prior art keywords
- image
- fusion
- nsst
- local contrast
- beta
- 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
- 239000011159 matrix material Substances 0.000 title claims abstract description 52
- 238000007500 overflow downdraw method Methods 0.000 title claims abstract description 17
- 230000004927 fusion Effects 0.000 claims abstract description 86
- 230000003595 spectral effect Effects 0.000 claims abstract description 49
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 22
- 230000000694 effects Effects 0.000 claims abstract description 15
- 230000009466 transformation Effects 0.000 claims abstract description 14
- 238000001228 spectrum Methods 0.000 claims abstract description 13
- 238000000034 method Methods 0.000 claims description 86
- 239000013598 vector Substances 0.000 claims description 12
- 238000007499 fusion processing Methods 0.000 claims description 9
- 230000003044 adaptive effect Effects 0.000 claims description 8
- 239000003086 colorant Substances 0.000 claims 1
- 239000002131 composite material Substances 0.000 claims 1
- 238000011156 evaluation Methods 0.000 abstract description 14
- 230000000007 visual effect Effects 0.000 abstract description 4
- 230000006870 function Effects 0.000 description 13
- 238000004458 analytical method Methods 0.000 description 8
- 238000000513 principal component analysis Methods 0.000 description 7
- 238000005457 optimization Methods 0.000 description 6
- 238000013135 deep learning Methods 0.000 description 4
- 238000012549 training Methods 0.000 description 4
- 230000009977 dual effect Effects 0.000 description 3
- 230000014759 maintenance of location Effects 0.000 description 3
- 238000001303 quality assessment method Methods 0.000 description 3
- 239000000654 additive Substances 0.000 description 2
- 230000000996 additive effect Effects 0.000 description 2
- 230000003190 augmentative effect Effects 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000013527 convolutional neural network Methods 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 238000003786 synthesis reaction Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012847 principal component analysis method Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 230000016776 visual perception Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- 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
- 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
- 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/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)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于NSST和β散度非负矩阵分解的遥感图像融合方法，属于多光谱和全色遥感图像融合技术领域，包括以下步骤：计算多光谱强度分量、光谱估计、NSST分解、低频系数融合、高频系数融合、NSST逆变换、图像重构。本发明采用上述一种基于NSST和β散度非负矩阵分解的遥感图像融合方法，解决了融合图像空间失真和光谱失真问题，在主观视觉效果和客观评价方面都达到了很好的效果，在提升空间分辨率的同时有效的保留了光谱信息。
Description
技术领域
本发明涉及多光谱和全色遥感图像融合技术领域，尤其是涉及一种基于NSST和β散度非负矩阵分解的遥感图像融合方法。
背景技术
多源遥感卫星可以提供众多具有不同空间分辨率、光谱分辨率和时间(时相)分辨率的遥感图像。高空间分辨率的全色图像(Panchromatic，PAN)，能够反映图像的整体空间结构信息，可详尽描述地物的细节特征。多光谱(Multispectral，MS)图像包含光谱信息，可以对地物进行识别、分类和解译，但是空间分辨率较低。将具有高空间分辨率的全色图像和光谱信息丰富的多光谱图像进行融合，能够获得具有较高空间分辨率的多光谱图像。这一过程又叫全色锐化，可以获得比单一类型图像更完整、更丰富的地表信息，从而改善后续处理效果。全色锐化广泛应用于土地利用规划、植被覆盖分析和地球资源调查等领域。
全色锐化方法可以分为成分替换法、多分辨率分析法、基于模型的方法三类。成分替换法通过线性或者非线性的变换将高光谱分辨率的多光谱图像变换到新的投影空间，分解为光谱成分和空间成分，并用全色图像去替换其空间成分，再经过逆变换获得融合图像。成分替换法主要包括IHS变换、主成分分析方法、GS方法、自适应GS(GSA)法、基于物理约束的波段相关空间细节(Band-dependent Spatial-detail with Physical Constrains,BDSD_PC)方法、基于部分置换的自适应分量替换(Partial Replacement-based AdaptiveComponent Substitution,PRACS)法等。成分替换法具有空间细节清晰、运行效率高、对误配准和混叠错误具有一定鲁棒性而被广泛使用，不过往往会带来光谱失真。
多分辨率分析方法是将全色图像多分辨率分解得到的空间细节注入到多光谱图像中，该方法通常通过线性分解方法获得，如小波变换、加性小波亮度比例(AdditiveWavelet Luminance Proportional,AWLP)方法等。与成分替换方法相比，多分辨率分析方法能更好地保持光谱特性，但容易产生空间结构失真。
基于模型的方法一般分为基于稀疏表示的方法和基于深度学习的方法。基于稀疏表示的方法首先从低空间分辨率数据中学习光谱字典，然后结合已知的高空间分辨率数据预测高空间分辨率和高光谱分辨率数据。例如，Li等提出了一种基于稀疏诱导先验信息的压缩感知方法，通过构建从多光谱图像中随机抽样的图像块字典来实现稀疏性。为了避免字典构建的成本，Zhu等提出了一种稀疏图像融合算法。Cheng等提出了一种融合小波变换和稀疏表示的融合框架。与多分辨分析法相比，这些方法具有超分辨能力和鲁棒性，能够获得更高的空间分辨率和光谱分辨率，且光谱失真较小。
近年来，人们对深度学习融合算法越来越感兴趣，例如Rao等提出了一种基于残差卷积神经网络的图像融合方法，直接学习输入与输出之间的残差，但浅层网络很难学习到深层次丰富的地物特征，数据量增多时还会出现欠拟合。Zhou利用深度学习技术增强对比度，合成全色图像，以在保持空间细节的同时减少光谱失真。Xiong等设计适用于全色锐化的损失函数和能够提取原始图像光谱和空间特征的四层卷积神经网络。Xiong等采用深度卷积神经网络学习全色图像和全色图像的光谱信息，并用光谱角控制光谱损失。Xu等提出基于模型的深度全色锐化方法。Xing等提出一种双协同融合模型。
虽然现有的融合算法在许多方面表现良好，但仍有一些方面需要改进。例如基于深度学习的方法往往需要非常大的训练集，而这些专业的遥感数据训练集通常很少。而且，不同的卫星有不同的数据类型，目前还很难做到不同卫星数据的同时训练。另外，网络训练需要大量的时间，使得网络调整的实时性较差。在基于稀疏表示的融合算法中，寻找最优的变换基比较困难。此外，稀疏表示有时会忽略图像的内在几何结构。
发明内容
本发明的目的是提供一种基于NSST和β散度非负矩阵分解的遥感图像融合方法，解决了融合图像空间失真和光谱失真问题，在主观视觉效果和客观评价方面都达到了很好的效果，在提升空间分辨率的同时有效的保留了光谱信息。
为实现上述目的，本发明提供了一种基于NSST和β散度非负矩阵分解的遥感图像融合方法，包括以下步骤：
S1、计算多光谱强度分量：使用基于加权局部对比度的自适应加权平均方法，融合多光谱图像各波段，生成强度分量I，加权局部对比度作为空间域中细节信息的评判指标，加权局部对比度高的像素被认为是权重更大的信息，在融合过程中赋予更大的权重，根据加权局部对比度设计自适应加权平均的系数ωi，公式如下：
其中n是MS图像的波段数，WSMLi表示MS图像第i个波段的加权局部对比度值；
S2、光谱估计：将I作为初始α，根据如下公式计算前景色F和背景色B：
其中i是第i颜色通道，αix和αiy值是光谱前景Fk、光谱背景Bk和α的水平和垂直导数，F和B包含光谱信息；
S3、NSST分解：对强度分量I和全色图像分别进行NSST分解，得到一个低频分量和多个高频分量，后续根据低频子带系数和各高频子带系数特点，实施不同的融合策略；
S4、低频系数融合：低频分量是原图像的近似，描述图像的基本结构，低频分量采用基于交替方向乘子法的β散度非负矩阵分解的融合规则；
S5、高频系数融合：NSST的不同尺度高频分量提供了多尺度信息，还包含边缘和纹理细节信息，边缘、纹理空间细节部分具有较高的局部对比度，是图像融合的目标，高频分量采用基于加权局部对比度的融合规则；
S6、NSST逆变换：对融合后的高、低频分量，进行NSST逆变换，得到融合图像，作为最终参与重构的α；
S7、图像重构：根据如下公式，通过组合α、F和B来重构得到最终的融合结果：
Ii＝αiFi+(1-αi)Bi (4)
其中Fi是第i个像素的前景颜色，Bi是第i个像素的背景颜色，Ii第i个像素的颜色，经过上面步骤之后，得到最终的融合结果。
优选的，所述步骤S2还包括图像抠图模型，具体为通过线性合成模型将输入图像区分为前景色F和背景色B，即第i个像素的颜色是对应的前景颜色和背景颜色的线性组合：
Ii＝αiFi+(1-αi)Bi (5)
其中Fi是第i个像素的前景色，Bi是第i个像素的背景色，α是F的不透明度，获取α是图像抠图的关键过程，根据图像抠图模型，在确定输入图像和α的同时，通过求解以下函数估计前景色F和背景色B：
其中i是第i颜色通道，αix和αiy值是光谱前景Fk、光谱背景Bk和α的水平和垂直导数。
优选的，所述步骤S5中加权局部对比度(WLCM)的具体步骤为NSST的不同尺度高频分量提供了多尺度信息，还包含边缘和纹理细节信息，在相同尺度下，边缘特征和纹理特征越明显，分量绝对值越大。但是，如果采用绝对值取最大作为高频分量的选择准则，会忽略邻近像素之间的相关性，也会给融合后的图像带来噪声，图像的边缘、纹理空间细节部分具有较高的局部对比度，是图像融合的目标，利用周围八个邻域的中值计算局部对比度，这样免高亮度噪声点被加权，从而被误判为细节信息，即减少了由高亮度孤立噪声引起的虚警，中心像素与周围像素的局部对比度为：
其中P0为局部区域中心像素的灰度值，Pmed为与中心像素相邻的8个邻域像素的中值，用下式计算：
Pmed＝median(Pi),i＝1,2,...,8. (9)
从上式中得出，如果中心区域A0为细节目标，即因此Cn＞P0，细节目标被增强；如果中心区域A0为背景，周围有细节信息或者都是背景，即/>因此Cn≤P0，背景被抑制；
若中心区域与周围邻域灰度差的均值越小，中心区域为空间细节的可能性越小，反之局部灰度差的均值越大，中心区域为空间细节的可能性越大，将中心像素与邻域的灰度差均值Mn作为局部对比度的权值，计算如下：
得到用加权的局部对比度公式：
WLCMn＝Cn*Mn. (10)。
优选的，所述步骤S4中基于交替方向乘子法的β散度非负矩阵分解为：
NMF问题的一般形式是：
两个矩阵之间的散度为元素散度之和：
β散度的表达式为：
引入新的变量W+和H+来应用非负性约束，约束条件为W＝W+和H＝H+，
重写为：
上面的符号表示由8个变量组成的增广拉格朗日函数，5个原始变量和3个对偶变量，从ADMM的角度来看，这是三部分优化：W、H和(X，W+，H+)，将优化目标分裂为X、W+和H+，对它们分别进行优化相当于共同优化它们：
在更新过程中，X的更新难度较大，且更新方法随参数β的取值情况而变化。
优选的，步骤S4中低频分量融合算法的具体为非负矩阵分解就是将非负矩阵分解为两个非负矩阵/>和/>相乘,使得X＝WH+ε，ε为背景噪声，此外，k<min{M，N}，原图像为真实图像在不同类型传感器中成像，加入一定背景噪声而获得的，即X＝WH+ε，将非负矩阵分解融合应用到MS和PAN图像融合过程中，在融合两种图像整体特征的基础上，保持多光谱图像的光谱特性，在低频分量的融合过程中，设置k＝1，首先，利用基于交替方向乘子法的β散度非负矩阵分解算法进行迭代，通过迭代使X和WH之间的重构误差达到最小，迭代完成后，得到唯一的特征基W，该矩阵包含参与融合的图像的总体特征，视为源图像的近似再现，使ε趋近于收敛，最后，将特征基W进行重置，还原为源图像的大小，即可获得效果很好的融合图像；
利用基于交替方向乘子法的β散度非负矩阵分解算法对低频分量AL和BL进行融合，具体的实现步骤为：
(1)将低频分量AL和BL分别按照行优先的方式整理成列向量的形式，得到列向量XA和XB，AL和BL的大小均为M×N，则列向量XA和XB的大小均为MN×l，具体如下所示：
(2)根据列向量XA和XB，构造原始数据矩阵X，其大小为MN×2；
(3)设定k＝1，NMF是有误差的分解，即X≈WH，定义某个目标函数，也称损失函数，目标函数衡量其逼近效果，目标函数选择KL(Kullback-Leibler)散度，随机生成初始迭代值W0和H0，设定最大迭代次数为2000，
W0＝rand(M,k),H0＝rand(k,N) (18)；
(4)设置好相关参数之后，利用基于交替方向乘子法的β散度非负矩阵分解算法，对原始数据矩阵X进行分解，迭代结束后，得到基矩阵W和权重系数矩阵H，W包含参与融合的低频分量AL和BL的整体特征，视为原图像的近似再现；
(5)将W进行重置变换，还原为M×N的矩阵L，L即为低频分量的融合图像。
优选的，步骤S5中高频分量融合算法具体为NSST的不同尺度高频分量提供了多尺度信息，还包含边缘和纹理细节信息，在相同尺度下，边缘特征和纹理特征越明显，分量绝对值越大，高频分量采用基于加权局部对比度的融合规则，具体如下：
MLCMD＝MLCMI(i,j)-MLCMP(i,j). (20)
其中，m、n分别是分解级数和方向数，代表像素点(i,j)处融合后的高频分量值，/>代表像素点(i,j)处强度分量I的高频系数值，/>代表代表像素点(i,j)处PAN图像的高频分量值，wI(i,j)为I的融合权重，wP(i,j)为全色图像的融合权重。
因此，本发明采用上述一种基于NSST和β散度非负矩阵分解的遥感图像融合方法，具备以下有益效果：
(1)本发明受图像抠图模型良好光谱保持性能的启发，将该模型引入到多光谱与全色图像的融合中，但是在遥感成像过程中，由于信噪比不同，多光谱和全色图像的特征并不完全相同。如果直接使用全色图像代替原α通道，会发生光谱畸变。所以，对传统的局部对比度进行改进，制定基于加权局部对比度的融合规则，根据融合规则，将多光谱图像各波段融合得到强度分量I。
(2)本发明分别对多光谱和全色图像进行NSST分解，可以得到一个低频分量和多个高频分量。在NSST分解的基础上，根据高低频分量信息的不同特点，设计针对其特点的融合规则。高频系数包含源图像丰富的边缘和纹理细节信息，采用加权局部对比度融合规则。低频分量是原图像的近似，描述图像的基本结构。本发明利用基于交替方向乘子法的β散度非负矩阵分解算法对低频分量进行融合。
(3)将全色图像与强度分量I融合后的图像作为新的α通道。根据图像抠图模型，根据前景色、背景色和α通道进行重构，最终可获得高空间分辨率和高光谱分辨率的融合图像。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明的融合流程图；
图2为中心像素及其八邻域示意图；
图3为三组MS和PAN图像；
图4为十种不同方法在第一组图像上的融合结果，其中(1)BT；(2)GSA；(3)GF；(4)IHS；(5)CT；(6)PCA；(7)PRACS；(8)NSST-SR；(9)BDSD-PC；(10)WT；(11)本发明的方法；(12)参考的多光谱图像；
图5为十种不同方法在第二组图像上的融合结果，其中(1)BT；(2)GSA；(3)GF；(4)IHS；(5)CT；(6)PCA；(7)PRACS；(8)NSST-SR；(9)BDSD-PC；(10)WT；(11)本发明的方法；(12)参考的多光谱图像；
图6为十种不同方法在第三组图像上的融合结果，其中(1)BT；(2)GSA；(3)GF；(4)IHS；(5)CT；(6)PCA；(7)PRACS；(8)NSST-SR；(9)BDSD-PC；(10)WT；(11)本发明的方法；(12)参考的多光谱图像。
具体实施方式
本发明提供了一种基于NSST和β散度非负矩阵分解的遥感图像融合方法，包括如下步骤：
(1)计算多光谱强度分量：
如果采用简单的平均融合规则，会丢失图像的细节信息。加权系数的准确选择决定了融合图像的质量。本发明使用基于加权局部对比度的自适应加权平均方法，融合MS图像各波段，生成强度分量I。加权局部对比度可以作为空间域中细节信息的评判指标。加权局部对比度高的像素被认为是权重更大的信息，如边缘或纹理等细节信息，在融合过程中赋予更大的权重。因此，本发明根据加权局部对比度设计自适应加权平均的系数ωi。
其中n是MS图像的波段数，WSMLi表示MS图像第i个波段的加权局部对比度值。
(2)光谱估计：
将I作为初始α，根据如下公式计算前景色F和背景色B。
其中i是第i颜色通道。αix和αiy值是光谱前景Fk、光谱背景Bk和α的水平和垂直导数。
F和B包含丰富的光谱信息，但不包含空间信息。后续步骤的主要目的是通过融合从PAN图像中获得空间细节信息。
(3)NSST分解：
对强度分量I和PAN图像分别进行NSST分解，得到相应的不同尺度和方向的分量。具体地说，可以得到一个低频分量和多个高频分量。本发明后续根据低频子带系数和各高频子带系数特点，实施不同的融合策略。
(4)低频系数融合：
低频分量是原图像的近似，只描述图像的基本结构，基本不包括边缘、轮廓等细节信息。本发明中低频分量采用基于交替方向乘子法的β散度非负矩阵分解的融合规则。
(5)高频系数融合：
NSST的不同尺度高频分量不仅提供了多尺度信息，还包含丰富的边缘和纹理细节信息。边缘、纹理等空间细节部分具有较高的局部对比度，是图像融合的目标。本发明中高频分量采用基于加权局部对比度的融合规则。
(6)NSST逆变换：
对融合后的高、低频分量，进行NSST逆变换，得到融合图像，作为最终参与重构的α。
(7)图像重构：
根据如下公式，通过组合α、F和B来重构得到最终的融合结果。
Ii＝αiFi+(1-αi)Bi (4)
经过上面步骤之后，便能得到最终的融合结果,本发明的流程图如图1。
补充说明：
1.图像抠图模型：
理论上，通过线性合成模型可以将输入图像区分为前景色F和背景色B，即第i个像素的颜色是对应的前景颜色和背景颜色的线性组合。
Ii＝αiFi+(1-αi)Bi (5)
其中Fi是第i个像素的前景色，Bi是第i个像素的背景色。α是F的不透明度。一般来说，获取α是图像抠图的关键过程。根据图像抠图模型，在确定输入图像和α的同时，可以通过求解以下函数估计前景色F和背景色B。
其中i是第i颜色通道。αix和αiy值是光谱前景Fk、光谱背景Bk和α的水平和垂直导数。
2.加权局部对比度：
NSST的不同尺度高频分量不仅提供了多尺度信息，还包含丰富的边缘和纹理细节信息。在相同尺度下，边缘特征和纹理特征越明显，分量绝对值越大。但是，如果采用绝对值取最大作为高频分量的选择准则，会忽略邻近像素之间的相关性，也会给融合后的图像带来噪声。中心像素及其八邻域如图2。
图像的边缘、纹理等空间细节部分具有较高的局部对比度，是图像融合的目标。传统的LCM算法用中心区域的像素最大灰度值的平方与周围八个邻域最大强度的比值来计算图像的局部对比度，这种方法易受高亮噪声影响，在运算过程中会引入虚警像素点，使得虚警率升高。本发明利用周围八个邻域的中值计算局部对比度。这样可以避免高亮度噪声点被加权，从而被误判为细节信息，即减少了由高亮度孤立噪声引起的虚警。因此，本发明定义中心像素与周围像素的局部对比度为:
其中P0为局部区域中心像素的灰度值。Pmed为与中心像素相邻的8个邻域像素的中值，可以用下式计算:
Pmed＝median(Pi),i＝1,2,...,8. (8)
从上式中可以得出，如果中心区域A0为细节目标，即因此Cn＞P0，细节目标被增强；如果中心区域A0为背景，周围可能有细节信息或者都是背景，即/>因此Cn≤P0，背景被抑制。
本发明初始利用邻域的中值计算局部对比度，避免了噪声点误判为空间细节。此外，本发明引入局部灰度差均值对局部对比度进行加权，计算加权的局部对比度，从而有效地实现对微弱空间细节的增强和对背景的抑制，较大程度的提高空间细节的显著性，提高细节信息检测率。若中心区域与周围邻域灰度差的均值越小，中心区域为空间细节的可能性越小，反之局部灰度差的均值越大，中心区域为空间细节的可能性越大。因此，本发明将中心像素与邻域的灰度差均值Mn作为局部对比度的权值，计算如下：
综上所述，可以得到用加权的局部对比度公式：
WLCMn＝Cn*Mn. (10)
3.基于交替方向乘子法的β散度非负矩阵分解：
NMF问题的一般形式是：
两个矩阵之间的散度被定义为元素散度之和：
β散度的表达式为：
W和H的非负性约束使W和H上的优化问题变得复杂。引入新的变量W+和H+来应用非负性约束，约束条件为W＝W+和H＝H+。
总之，可以重写为：
上面的符号表示由8个变量组成的增广拉格朗日函数，5个原始变量和3个对偶变量。从ADMM的角度来看，这是三部分优化：W、H和(X，W+，H+)。这是因为将优化目标分裂为X、W+和H+，因此对它们分别进行优化相当于共同优化它们：
在更新过程中，X的更新难度较大，且更新方法随参数β的取值情况而变化。
4.低频分量融合算法：
非负矩阵分解就是将非负矩阵分解为两个更小的非负矩阵/>和/>相乘，使得X＝WH+ε，ε为背景噪声。此外，k远小于M和N，即k<min{M，N}。原图像通常可以视为真实图像在不同类型传感器中成像，加入一定背景噪声而获得的，即X＝WH+ε。将非负矩阵分解融合应用到MS和PAN图像融合过程中，可以在融合两种图像整体特征的基础上，尽可能多地保持多光谱图像的光谱特性，从而达到遥感图像融合的目的。
本发明在低频分量的融合过程中，设置k＝1。首先，利用基于交替方向乘子法的β散度非负矩阵分解算法进行迭代。迭代求解实际上是一个优化过程。通过将迭代使X和WH之间的重构误差达到最小，可以有效抑制背景噪声。迭代完成后，可以得到唯一的特征基W，该矩阵包含参与融合的图像的总体特征，视为源图像的近似再现。可以使ε趋近于收敛，并且有效抑制背景噪声。最后，将特征基W进行重置，还原为源图像的大小，即可获得效果很好的融合图像。
本发明利用基于交替方向乘子法的β散度非负矩阵分解算法对低频分量AL和BL进行融合。具体的实现步骤为:
(1)将低频分量AL和BL分别按照行优先的方式整理成列向量的形式，得到列向量XA和XB。AL和BL的大小均为M×N，则列向量XA和XB的大小均为MN×l。具体如下所示。
(2)根据列向量XA和XB，构造原始数据矩阵X，其大小为MN×2。
(3)设定k＝1。NMF是有误差的分解，即X≈WH。为了获得一个近似的分解，使得X和WH之间的重构误差最小，必须定义某个目标函数，也称损失函数。目标函数可以衡量其逼近效果。本方法目标函数选择KL(Kullback-Leibler)散度。随机生成初始迭代值W0和H0，本发明设定最大迭代次数为2000。
W0＝rand(M,k),H0＝rand(k,N) (18)
(4)设置好相关参数之后，利用基于交替方向乘子法的β散度非负矩阵分解算法，对原始数据矩阵X进行分解。迭代结束后，得到基矩阵W和权重系数矩阵H。W包含参与融合的低频分量AL和BL的整体特征，可以视为原图像的近似再现。
(5)将W进行重置变换，还原为M×N的矩阵L。L即为低频分量的融合图像。
4.高频分量融合算法：
NSST的不同尺度高频分量不仅提供了多尺度信息，还包含丰富的边缘和纹理细节信息。在相同尺度下，边缘特征和纹理特征越明显，分量绝对值越大。但是，如果采用绝对值取最大作为高频分量的选择准则，会忽略邻近像素之间的相关性，也会给融合后的图像带来噪声。本发明中高频分量采用基于加权局部对比度的融合规则，具体如下：
MLCMD＝MLCMI(i,j)-MLCMP(i,j). (20)
其中，m、n分别是分解级数和方向数，代表像素点(i,j)处融合后的高频分量值，/>代表像素点(i,j)处强度分I的高频系数值，/>代表代表像素点(i,j)处PAN图像的高频分量值。wI(i,j)为I的融合权重，wP(i,j)为PAN图像的融合权重。
实施例
为了说明本发明的有效性，特进行如下实验论证：
本发明使用了一个包含27组图像的数据集进行实验，这些图像是由LANDSAT 7ETM+拍摄的，在6个波段(红、绿、蓝、近红外、中红外和短波红外)工作。多光谱图像的空间分辨率为30m，全色图像的空间分辨率为15m。因此，由于数据集中没有高分辨率的MS图像作为参考图像，所以我们先将原始的多光谱图像进行上采样，得到像素大小为400×400的多光谱图像。接着，将像素大小为400×400的多光谱图像和全色图像进行下采样，得到像素大小为200×200的多光谱图像和全色图像作为实验图像。随机选择三对不同场景的图像进行对照实验。最后，将原始MS图像作为参考图像，与各方法融合后的图像进行比较。图3展示了三组MS和PAN图像，随后将用于实验分析。
本发明将十种具有代表性的融合方法：(1)BT(Brovey Transform-basedmethod)；(2)GSA(Gram Schmidt Adaptive-based method)；(3)GF(Guided Filter-basedmethod)；(4)IHS(Intensity-Hue-Saturation-based method)；(5)CT(CurveletTransform-based method)；(6)PCA(Principal Component Analysis-based method)；(7)PRACS(Partial Replacement Adaptive Component Substitution-based method)；(8)NSST-SR(NSST and Sparse Representation-based method)；(9)BDSD-PC(Banddependent Spatial-detail with Physical Constrains-basedmethod)；(10)WT(WaveletTransform-based method)和本发明的融合方法进行比较。
通常可以通过主观评价和客观评价来衡量一种遥感图像融合方法的性能。在进行主观评价时，通常会考虑目标的清晰度以及融合图像与原始多光谱图像光谱的接近程度。然而，仅凭主观评价难以对融合质量进行准确比较。为了对图像融合方法进行定量评价，采用若干指标来评价不同融合方法的性能。在实验中，使用六个众所周知的客观评价指标，并详细介绍如下。
(1)相关系数(Correlation Coefficient，CC)计算参考图像和融合结果之间的相关性。该值越大，表示融合结果越接近参考图像，理想值是1。
(2)光谱角制图(SpectralAngle Mapper，SAM)反映融合图像和参考图像之间的光谱失真。SAM值越小，融合图像中的光谱失真越小。
(3)光谱信息散度(Spectral Information Divergence，SID)评测光谱之间的差异，理想值为0。
(4)无参考质量评价指标(Quality with No Reference，QNR)可以在没有参考图像的情况下评估融合图像的质量，由光谱失真指数Dλ、空间失真指数DS和全局QNR值三部分组成。对于全局QNR，值越大则融合效果越好，理想值为1。
(5)Dλ是QNR的亚度量，可以测量光谱失真。其值越小，融合效果越好，理想值为0。
(6)DS是QNR的亚度量，可以测量空间畸变。其值越小，融合效果越好，理想值为0。
图4、5、6中的(1)-(11)分别给出第一组图像、第二组图像、第三组图像分别经(1)BT；(2)GSA；(3)GF；(4)IHS；(5)CT；(6)PCA；(7)PRACS；(8)NSST-SR；(9)BDSD-PC；(10)WT方法与本发明的方法得到的融合结果，并与参考的多光谱图像(12)作对比。为了更加直观地比较融合结果之间的差异，还对融合结果细节进行局部放大，放大图像置于图像中右下角。表1、2、3分别列出三组图像融合结果的六种包含光谱和空间质量评价的客观指标数值，对于所有的融合质量评价标准，最好的结果以红色粗体显示。
表1第一组图像融合结果的客观评价
表2第二组图像融合结果的客观评价
表3第三组图像融合结果的客观评价
从图4，5，6中可以看出，BT、GSA、GF、IHS方法的融合图像空间细节比较清晰，比较完整的保留了全色图像的空间细节信息。但是在整体区域出现光谱失真，在局部放大区域比较明显。PCA方法的融合图像空间细节和光谱均发生了严重失真。WT方法的融合图像出现了虚影，空间细节失真较为严重。其余方法保持了较好的光谱特性，但是局部放大部分空间细节比较模糊。从主观视觉角度分析，本发明的空间细节更加清晰，在改善空间细节的同时具有很好的光谱保持特性。
从表1、2、3中可以看出，本发明在六个评价指标中的CC、SID、Ds、QNR上表现最好。尤其是QNR值，在所有方法中遥遥领先，比较接近最优值1。Dλ和SAM值也都较小，和最优方法相比差别很小。从客观评价指标角度分析，本发明具有良好的空间细节保持特性和光谱保持特性，整体效果较好。
综上所述，本发明在视觉感知和客观评价方面都能达到较好的效果。能够在保留多光谱图像更多光谱信息的同时，从全色图像中获得更多的空间细节，验证了本发明的有效性。
因此，本发明采用上述一种基于NSST和β散度非负矩阵分解的遥感图像融合方法，解决了融合图像空间失真和光谱失真问题，在主观视觉效果和客观评价方面都达到了很好的效果，在提升空间分辨率的同时有效的保留了光谱信息。
最后应说明的是：以上实施例仅用以说明本发明的技术方案而非对其进行限制，尽管参照较佳实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对本发明的技术方案进行修改或者等同替换，而这些修改或者等同替换亦不能使修改后的技术方案脱离本发明技术方案的精神和范围。
Claims (1)
1.一种基于NSST和β散度非负矩阵分解的遥感图像融合方法，其特征在于，包括以下步骤：
S1、计算多光谱强度分量：使用基于加权局部对比度的自适应加权平均方法，融合多光谱图像各波段，生成强度分量I，加权局部对比度作为空间域中细节信息的评判指标，加权局部对比度高的像素被认为是权重更大的信息，在融合过程中赋予更大的权重，根据加权局部对比度设计自适应加权平均的系数，公式如下：
（1）；
（2）；
其中n是MS图像的波段数，WSML i 表示MS图像第i个波段的加权局部对比度值；
S2、光谱估计：将I作为初始α，根据如下公式计算前景色F和背景色B：
（3）；
其中是第/>颜色通道，/>、/>、/>、/>、/>和/>值是光谱前景/>、光谱背景/>和/>的水平和垂直导数，F和B包含光谱信息；
S3、NSST分解：对强度分量I和全色图像分别进行NSST分解，得到一个低频分量和多个高频分量，后续根据低频子带系数和各高频子带系数特点，实施不同的融合策略；
S4、低频系数融合：低频分量是原图像的近似，描述图像的基本结构，低频分量采用基于交替方向乘子法的β散度非负矩阵分解的融合规则；
S5、高频系数融合：NSST的不同尺度高频分量提供了多尺度信息，还包含边缘和纹理细节信息，边缘、纹理空间细节部分具有较高的局部对比度，是图像融合的目标，高频分量采用基于加权局部对比度的融合规则；
S6、NSST逆变换：对融合后的高、低频分量，进行NSST逆变换，得到融合图像，作为最终参与重构的；
S7、图像重构：根据如下公式，通过组合、F和B来重构得到最终的融合结果：
（4）；
其中是第/>个像素的前景颜色，/>是第/>个像素的背景颜色，/>第/>个像素的颜色，经过上面步骤之后，得到最终的融合结果；
所述步骤S5中加权局部对比度的具体步骤为NSST的不同尺度高频分量提供了多尺度信息，还包含边缘和纹理细节信息，在相同尺度下，边缘特征和纹理特征越明显，分量绝对值越大，如果采用绝对值取最大作为高频分量的选择准则，会忽略邻近像素之间的相关性，也会给融合后的图像带来噪声，图像的边缘、纹理空间细节部分具有较高的局部对比度，是图像融合的目标，利用周围八个邻域的中值计算局部对比度，这样免高亮度噪声点被加权，从而被误判为细节信息，即减少了由高亮度孤立噪声引起的虚警，中心像素与周围像素的局部对比度为：
（5）；
其中P 0 为局部区域中心像素的灰度值，P med 为与中心像素相邻的8个邻域像素的中值，用下式计算：
（6）；
从上式中得出，如果中心区域A 0 为细节目标，即，因此/>，细节目标被增强；如果中心区域A 0 为背景，周围有细节信息或者都是背景，即/>，因此/>，背景被抑制；
若中心区域与周围邻域灰度差的均值越小，中心区域为空间细节的可能性越小，反之局部灰度差的均值越大，中心区域为空间细节的可能性越大，将中心像素与邻域的灰度差均值M n 作为局部对比度的权值，计算如下：
（7）；
得到用加权的局部对比度公式：
（8）；
步骤S4中低频分量融合算法的具体为非负矩阵分解就是将非负矩阵，分解为两个非负矩阵/>和/>相乘,使得/>，/>为背景噪声，此外，k<min{M，N}，原图像为真实图像在不同类型传感器中成像，加入一定背景噪声而获得的，即，将非负矩阵分解融合应用到MS和高空间分辨率的全色图像PAN融合过程中，在融合两种图像整体特征的基础上，保持多光谱图像的光谱特性，在低频分量的融合过程中，设置k=1，首先，利用基于交替方向乘子法的β散度非负矩阵分解算法进行迭代，通过迭代使X和WH之间的重构误差达到最小，迭代完成后，得到唯一的特征基W，该矩阵包含参与融合的图像的总体特征，视为源图像的近似再现，使/>趋近于收敛，最后，将特征基W进行重置，还原为源图像的大小，即可获得效果很好的融合图像；
利用基于交替方向乘子法的β散度非负矩阵分解算法对低频分量AL和BL进行融合，具体的实现步骤为：
（1）将低频分量AL和BL分别按照行优先的方式整理成列向量的形式，得到列向量XA和XB，AL和BL的大小均为M×N，则列向量XA和XB的大小均为MN×l，具体如下所示：
（9）；
（2）根据列向量XA和XB，构造原始数据矩阵X，其大小为MN×2；
（10）；
（3）设定k=1，NMF是非负矩阵，NMF是有误差的分解，即X≈WH，定义某个目标函数，也称损失函数，目标函数衡量其逼近效果，目标函数选择KL（Kullback-Leibler）散度，随机生成初始迭代值W 0 和H 0 ，设定最大迭代次数为2000，
（11）；
（4）设置好相关参数之后，利用基于交替方向乘子法的β散度非负矩阵分解算法，对原始数据矩阵X进行分解，迭代结束后，得到基矩阵W和权重系数矩阵H，W包含参与融合的低频分量AL和BL的整体特征，视为原图像的近似再现；
（5）将W进行重置变换，还原为M×N的矩阵L，L即为低频分量的融合图像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210657229.7A CN114926452B (zh) | 2022-06-10 | 2022-06-10 | 一种基于NSST和β散度非负矩阵分解的遥感图像融合方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210657229.7A CN114926452B (zh) | 2022-06-10 | 2022-06-10 | 一种基于NSST和β散度非负矩阵分解的遥感图像融合方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114926452A CN114926452A (zh) | 2022-08-19 |
| CN114926452B true CN114926452B (zh) | 2024-04-02 |
Family
ID=82814974
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210657229.7A Active CN114926452B (zh) | 2022-06-10 | 2022-06-10 | 一种基于NSST和β散度非负矩阵分解的遥感图像融合方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114926452B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115861141B (zh) * | 2022-12-02 | 2023-06-16 | 北京领云时代科技有限公司 | 基于pcnn神经网络的无人机获取图像处理系统及方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107169946A (zh) * | 2017-04-26 | 2017-09-15 | 西北工业大学 | 基于非负稀疏矩阵与超球面彩色变换的图像融合方法 |
| CN114240990A (zh) * | 2021-12-07 | 2022-03-25 | 电子科技大学 | Sar图像点目标分割方法 |
-
2022
- 2022-06-10 CN CN202210657229.7A patent/CN114926452B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107169946A (zh) * | 2017-04-26 | 2017-09-15 | 西北工业大学 | 基于非负稀疏矩阵与超球面彩色变换的图像融合方法 |
| CN114240990A (zh) * | 2021-12-07 | 2022-03-25 | 电子科技大学 | Sar图像点目标分割方法 |
Non-Patent Citations (4)
| Title |
|---|
| A Closed-Form Solution to Natural Image Matting;Anat Levin 等;IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE;20080229;第30卷(第2期);第1-15页 * |
| ALTERNATING DIRECTION METHOD OF MULTIPLIERS FOR NON-NEGATIVE MATRIX FACTORIZATION WITH THE BETA-DIVERGENCE;Dennis L. Sun 等;2014 IEEE International Conference on Acoustic, Speech and Signal Processing (ICASSP);20141231;第1-5页 * |
| 一种NSST与稀疏表示相结合的遥感图像融合算法;金益如;杨学志;董张玉;郑鑫;李国强;;地理与地理信息科学;20160315(第02期);全文 * |
| 结合视觉显著性与Dual-PCNN的红外与可见光图像融合;侯瑞超;周冬明;聂仁灿;刘栋;郭晓鹏;;计算机科学;20180615(第S1期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114926452A (zh) | 2022-08-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108830796B (zh) | 基于谱空结合和梯度域损失的高光谱图像超分辨重构方法 | |
| CN111047516B (zh) | 图像处理方法、装置、计算机设备和存储介质 | |
| Song et al. | Unsupervised denoising for satellite imagery using wavelet directional CycleGAN | |
| EP2230855B1 (en) | Synthesizing virtual images from texture and depth images | |
| US20080267525A1 (en) | Soft Edge Smoothness Prior and Application on Alpha Channel Super Resolution | |
| CN110189286B (zh) | 一种基于ResNet的红外与可见光图像融合方法 | |
| CN114897882B (zh) | 一种基于加权平均曲率滤波分解的遥感图像融合方法 | |
| CN109933639B (zh) | 一种面向图层叠加的多光谱图像与全色图像自适应融合方法 | |
| CN114926452B (zh) | 一种基于NSST和β散度非负矩阵分解的遥感图像融合方法 | |
| CN115760814A (zh) | 一种基于双耦合深度神经网络的遥感图像融合方法及系统 | |
| CN115861083B (zh) | 一种多尺度和全局特征的高光谱与多光谱遥感融合方法 | |
| CN116645569A (zh) | 一种基于生成对抗网络的红外图像彩色化方法和系统 | |
| CN110830043B (zh) | 一种基于混合加权全变分和非局部低秩的图像压缩感知重构方法 | |
| CN114897694A (zh) | 基于混合注意力和双层监督的图像超分辨率重建方法 | |
| CN112270646A (zh) | 基于残差密集跳跃网络的超分辨增强方法 | |
| CN116739899A (zh) | 基于saugan网络的图像超分辨率重建方法 | |
| CN116309136A (zh) | 一种基于sar先验知识引导的遥感影像云区重建方法 | |
| CN114638749A (zh) | 低光照图像增强模型、方法、电子设备及存储介质 | |
| CN112200752B (zh) | 一种基于er网络多帧图像去模糊系统及其方法 | |
| Heimann et al. | Frequency-selective mesh-to-mesh resampling for color upsampling of point clouds | |
| Ponomaryov et al. | Fuzzy color video filtering technique for sequences corrupted by additive Gaussian noise | |
| CN114897757A (zh) | 一种基于nsst和参数自适应pcnn的遥感图像融合方法 | |
| Song et al. | Unsupervised denoising for satellite imagery using wavelet subband cyclegan | |
| Liang et al. | Multi-scale and multi-patch transformer for sandstorm image enhancement | |
| CN113724139B (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 |