IHS变换 python ihs变换 为什么用分量I

前言

在遥感图像融合领域，说到已经投入到商业运用，并且其变种还十分多的融合算法，恐怕IHS算法要算其中之一了。对于初学者，比如我，也经常在论文中看到IHS及其变种的身影。故在此就目前自己对IHS算法的了解做个笔记。

IHS算法简介

IHS即是Intensity–Hue–Saturation的缩写，同我们熟知的RGB一样，也是一种颜色空间。3者的意思为，强度，色调，饱和度。通常IHS算法只对强度分量I进行操作，而无关其他2个分量。确切地说，进行IHS变换之后，强度分量I只包含图像的结构信息而去除了其他的信息，剔除掉的信息主要是颜色信息。而IHS所属的成分替代方法假设Pan图所包含的结构同强度分量I是一样的，所以实践中能够使用Pan图直接替换强度分量I得到融合结果。
那么如何计算强度分量I呢？由于IHS算法假设(assumption)强度分量I是多光谱图像(MS)各波段的线性组合，所以我们可以使用低分辨MS(LRM)图像各波段的线性组合得到强度分量I。虽然根据光谱响应函数来看，即使Pan图与多光谱各波段存在某种关系，也绝对不是线性的。换句话说，假设强度分量I是多光谱各波段的线性组合是不精确的，但是由于其在建模和求解上更易实现，而且效果也不错，所以仍不失为一种解决问题的方法。其特点是计算效率很快，融合结果的空间分辨率好，但是光谱失真(color distortion)严重。部分原因可在【2006】《A New Intensity Hue-Saturation Fusion Approach to Image Fusion With a Tradeoff Parameter》一文中找到，Choi指出：IHS算法的光谱失真严重主要是因为Pan和强度分量I之间差异过大，所以通过$min_{I_{new}}\{|Pan-I_{new}|^2+|I_{new}-I|^2\}$获得一个新的强度分量。其结果就是在光谱保持和空间细节注入方面取得权衡。当$I_{new}$更加逼近Pan时，融合结果的空间分辨率更高，当然光谱失真就愈严重。而当$I_{new}$更逼近I时，则融合结果的光谱保持更好，但是空间分辨率便不理想。需要注意的是，原始的IHS算法只适用于RGB3个通道的图像，后来有很多IHS的扩展算法将其扩展到4通道乃至n通道。最一般的确定n通道MS图像对应的强度分量的I的公式如下。$I=∑_{i=1}^{N}α_{i}M_{i} (1)$
其中$α_{i}$表示各波段的系数，$M_{i}$表示上采样的LRM，N表示MS的通道总数。通常我们都取$α_{i}=\frac{1}{N}$ 。当然也有更加精确的确定I的公式，比如choi在【2006】年发表的《A new intensity-hue-saturation fusion approach to image fusion with a tradeoff parameter》中指出强度分量I等于$I=\frac{1}{3}(0.3*R+0.75*G+0.25*B+1.7*NIR)$；而Tu在[2004]《A fast intensity-hue-saturation fusion technique with spectral adjustment for IKONOS imagery》[FIHS模型]一文中指出IKONOS中的强度分量I为$I=\frac{1}{3}(R+0.75*G+0.25*B+NIR)$；对于QuickBird卫星图像融合，Joshi建议强度分量I应当等于$I=0.2308*R+0.2315*G+0.1139*B+0.4239*NIR$《A model-based approach to multiresolution fusion in remotely sensed images》【2006】。当然也有使用自适应参数算法求得$I=∑_{i=1}^{N}α_{i}M_{i}$中的$α_{i}$，比如在【2010】年发表的《An Adaptive IHS Pan-Sharpening Method》基于Pan尽可能逼近强度分量可减少光谱失真的假设求得$α_{i}$，即$P≈∑_{i=1}^{N}α_{i}M_{i}$。然后利用这些$α_{i}$作为权重系数求出强度分量I，最后使用由Pan图决定的权重矩阵式(2)改进FIHS模型式(3)进行细节注入，但是光谱失真仍然存在。 $h(x)=exp(-\frac{λ}{|▽P|^4+ε}) (2)$
其中$λ=10^{-9}$， $ε=10^{-10}$。
$F_{i} = M_{i} + h(x)(P − I) (3)$
在【2014】年《An Improved Adaptive Intensity–Hue–Saturation Method for the Fusion of Remote Sensing Images》一文中，由于使用Pan图的权重矩阵引入的权重过大，通过添加MS的权重矩阵式(4)对FIHS模型进行改进得到式(5)以决定细节注入。虽然改善了光谱失真，但是导致细节丰富的区域的结果过于平滑，图像细节较AIHS算法略平滑。
$W_{M_{i}}=exp(-\frac{λ}{|▽M_{i}|^4+ε}) (4)$
$F_{i} = M_{i} + W_{i}(P − I) (5)$
这里通过权衡Pan的权重矩阵和多光谱权重矩阵的影响，得到一个折中的权重矩阵以控制细节的注入。$W_{i}$的计算公式如下。
$W_{i} =\frac{M_{i}}{ \frac{1}{N}∑_{i=1}^{N}M_{i}} (βW_{P}+(1-β)W_{M_{i}}) (6)$
其中$W_{P}$就是前面提到的式(2)，只不过在此为了命名上的统一称为$W_{P}$。

IHS结合变分模型

近年来对于IHS结合变分模型的算法也取得了良好的效果，在一定程度上减少了光谱失真，如《Panchromatic and Multi-spectral Image Fusion Using IHS and Variational Models》【2012】一文使用变分法求得替代成分f，不再使用Pan直接替代求得强度分量，而是将f，H，S应用在IHS逆变换中，实验证明该方法在光谱保持上取得了一定的进步。其主要思想就是向f中注入Pan图像的细节，并且通过逼近I来保持光谱信息。在《JOINT IHS AND VARIATIONAL METHODS FOR PAN-SHARPENING OF VERY HIGH RESOLUTION IMAGERY》【2013】Zhou提出由3项组成的能量函数来求得替代成分的最优解。他认为应该从Pan图提取空间信息注入到MS中以添加细节，通过使替代成分f的低分辨版本接近于强度成分I保持光谱特征，第三项用于保证解的平滑性。在【2016】年Zhou又发表《EXTENDED GIHS FUSION FOR PAN-SHARPENING BASED ON IMAGE MODEL》实现了IHS和变分结合的模型。另外Ghahremani于2016年提出一种新的基于变分的模型NIHS(非线性的IHS)《Nonlinear IHS: A Promising Method for Pan-Sharpening》。通过引入局部合成模型使用非线性的方法近似强度分量，引入全局模型使高空间分辨率下的强度分量与其退化(即下采样的强度分量)的强度分量协调一致。具体来说，通过将Pan图及其下采样版本和LRM上采样版本及LRM分别分成对应大小的patch，然后将Pan图及其下采样版本的patch看做LRM上采样版本及LRM各波段对应位置的patch的线性组合实现非线性合成方法。值得注意的是该局部合成方法同样利用了Pan图应当尽可能地逼近强度分量可以减轻光谱失真的假设。

其中i表示patch的序号，$s^{(i)},S^{(i)}$分别表示LRM的强度分量及其上采样版本的第i个patch，$Y^{(i)} =( y_{1}^{(i)} , ....,y_{k}^{(i)} )^T$和$Y_{up}^{(i)} =( y_{1,up}^{(i)} , ....,y_{k,up}^{(i)} )^T$分别表示LRM及其上采样版本在波段方向的第i个patch，k表示波段序号，个人认为这里的L表示波段总数。

其中$X^{(i)},x^{(i)}$分别表示Pan及其下采样版本的第i个patch。

然后全局模型利用局部合成方法求出的中间结果得到满足全局重构的最终结果，即新的强度分量。全局模型使用$I_{up} =argmin_{I_{up}}\{||I-MI_{up}||_{2}^2+η||I_{up}-I_{0,up}||_{2}^2\} (10)$

其中$I_{up}$为最终得到的强度分量，而$I和I_{0,up}$分别为通过局部合成方法得到的强度分量及其上采样版本。

实验结果表明该算法在保持光谱信息方面表现优秀，空间分辨率也尚可，但是在边缘处略有模糊会出现块状阴影。就普适性和客观评价结果来看，该算法的效果相当不错。

IHS和小波结合的模型

借助于小波变换能够很好地保持原始多光谱图像的光谱信息，而IHS算法通常能够平滑地融合色彩和空间特征的优点。2004年，zhang在《An IHS and wavelet integrated approach to improve pan-sharpening visual quality of natural colour IKONOS and QuickBird images》一文中提出使用IHS结合小波变换模型进行多光谱图像融合。通过对强度分量I和直方图匹配之后的Pan图分别进行一层的小波分解，使用强度分量I和Pan图的低频子图的线性组合替换Pan的低频子图以保持光谱信息，然后对该中间结果实施小波重构得到一幅新的图像当做Pan图，最后使用新的Pan图代替强度分量I。这里介绍一下如何计算新的低频子图，$LL_{new} =w_{1}*I_{LL}+w_{2}*P_{LL}$，其中权重系数使用下面的公式计算。

这里$\hat{a}$,$\hat{b}$分别表示a,b的均值，N表示近似图像的像素点总数，即$I_{LL}$或$P_{LL}$(因为2者一样大)的像素点总数。

算法流程图如下。

下面就来介绍一下IHS融合算法，并使用matlab实现IHS算法。

IHS算法步骤

首先我们按照最传统的IHS算法的讲解方式来讲解IHS算法，按照该思路，IHS算法可以分为4步：1）对上采样的MS做IHS正变换；2）对Pan图进行直方图匹配得到$Pan_{new}$；3）使用$Pan_{new}$替代强度分量I；4）对替代之后的结果进行IHS逆变换，最终得到融合的结果。
以上就是IHS算法的全部内容了，需要思考的是，为什么要对Pan进行直方图匹配？An Improved Adaptive Intensity–Hue–Saturation Method for the Fusion of Remote Sensing Images【2014】一文中指出，直方图匹配主要是为了消除大气、光照和传感器不同的影响。
接下来在数学的角度来看一看，IHS是如何实现的，依据此分析便可以实现IHS的融合算法了。

IHS算法公式形式

1）IHS变换的公式就如下所示

2）直方图匹配，你既可以使用matlab中imhist，histeq来进行，也可以使用下面的直方图匹配公式，其中σ表示方差，μ表示均值。

3）然后我们使用上面的Pan_new来替代IHS中的I

4）接着使用逆变换公式来得到融合结果

Fast IHS算法

如果没有看懂上面的IHS算法也没有关系，实际上IHS算法可以更简洁地实现，这得益于Tu在2004年的研究工作。使得像我这样的渣渣对IHS算法的理解和实现都容易了许多。像上面的步骤实现IHS还是有点繁琐，当我们使用公式求得强度分量，再将Pan图直方图匹配到强度分量，就可以直接使用下面的公式实现IHS算法。至于严苛的数学证明这里就不再赘述。

可以看到，除了Pan依旧要做直方图匹配，IHS的逆变换可以绕过，该方法被称为Fast IHS算法。由Tu在2004年提出，当然该模型可以很容易扩展到n通道，公式如下：$M_{i}^{f}=M_{i}+(P-I)$

matlab代码

鉴于显示便于观察融合效果的缘故，我们使用3通道的图像作为演示IHS融合算法的实例。
主函数代码如下：

LMS = im2double(imread('LR.tif'));
MS = imresize(LMS,4);
Pan = im2double( imread('pan.tif') );
fused = IHS(MS,Pan);
figure(1);imshow(fused);title('IHS融合结果');

IHS算法代码如下：注意MS各波段前的系数α按照最一般的设置$\frac{1}{N}$来取的，但也有通过改变该系数形成的方式来改善IHS算法的。如：An Adaptive IHS Pan-Sharpening Method【2010】，文中通过求解自适应系数α来获取更加精确的I，使得Pan与I的差异减小，文中指出之所以IHS算法光谱失真严重，主要就在于Pan与I之间的差异太大，导致注入过多的空间信息。

function [fused] = IHS(MS,Pan)
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % 输入参数：
% %        MS：上采样之后的多光谱图像
% %        Pan ： 全色图
% % 输出参数：
% %        fused: IHS融合结果
% % 
% %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
   I = zeros(size(MS));
   for n=1:3
        I = I + MS(:,:,n)/n ; 
    end
    
    P = (Pan - mean(Pan(:))) * std(I(:))/std(Pan(:)) + mean(I(:)); % 直方图匹配
    fused = zeros(size(MS));
    
    for n=1:3
        fused(:,:,n)= MS(:,:,n) + P - I; 
    end
end

融合结果

今天就先写到这，以后还有新的了解就再添上。