CN105354799A - 一种盲目反卷积红外光谱图超分辨率复原方法 - Google Patents

Abstract

本发明公开了一种盲目反卷积红外光谱图超分辨率复原方法，该方法包括步骤：对红外光谱数据进行归一化处理；计算归一化后的光谱强度的一阶导数；根据光谱一阶导数估计出光谱的噪声强度，从而确定正则化参数的大小；建立红外光谱模型；利用Huber-Markov场、光谱反卷积和交替最小值最优迭代算法求解高分辨率的红外光谱。本发明在光谱平坦区域加强噪声抑制，在光谱陡峭区域加强光谱结构保存，对于所用的傅里叶红外光谱，空域迭代盲目反卷积算法效果良好，在信噪比高于30dB时，光谱分辨率提高约30%；本发明提出的方法有效利用了Huber-Markov自适应约束的先验知识，促进了红外光谱在目标检测识别及农业、生物、医学、气象、空间技术、石油化工领域的应用，并提供技术支撑。

Description

一种盲目反卷积红外光谱图超分辨率复原方法

技术领域

本发明涉及红外光谱数据处理技术领域，尤其涉及一种盲目反卷积红外光谱图超分辨率复原方法。

背景技术

近几十年，人们对中波红外光谱技术的研究相当活跃，使得红外光谱在目标跟踪、化学计量分析、材料分析和生学医学等领域的得到了广泛的应用。然而，红外光谱数据在光谱仪测量过程中，经常会受到谱线的碰撞展宽、自然展宽、电路和探测器低通特性等因素的影响，导致光谱数据出现谱带重叠和强噪声干扰现象。这些干扰给红外光谱应用带来了困难，甚至会导致严重的物质鉴定和目标识别误差，极大限制了红外光谱在实际生活中的应用。因此，对光谱数据进行预处理、提高光谱的分辨率就显得非常重要。

目前，复原退化光谱的方法有很多，这些方法主要可分为两大类：光谱反卷积和曲线拟合技术。曲线拟合技术认为光谱曲线是通过高斯函数、洛伦兹函数或它们的组合而构成，即实测光谱能通过这些函数来拟合得到。这种方法能够很好的解决谱带重叠的问题，即能拟合出重叠的波峰，但是对光谱噪声极其敏感。光谱反卷积技术主要是通过不同的约束来复原退化的光谱，这些方法有傅里叶自反卷积(FourierSelf-deconvolution，FSD)，最大香农熵反卷积、Tikhonov正则化、同态滤波等方法。

发明内容

本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种盲目反卷积红外光谱图超分辨率复原方法，解决了由傅里叶红外光谱仪测得的红外光谱数据经常会受到随机噪声和仪器误差等的影响而导致光谱分辨率降低的问题。

本发明解决其技术问题所采用的技术方案是：一种盲目反卷积红外光谱图超分辨率复原方法，包括以下步骤：

1)采集待处理的红外光谱，对待处理的红外光谱数据作归一化处理，其中i表示波数点，n表示波数的点数；

2)对归一化后光谱强度f_i进行一阶求导，f_i′＝(f_i+1-f_i)/2；

3)根据光谱一阶导数强度值f′_i，使用绝对中值估计器初步计算光谱噪声水平，并由此确定正则化参数的数值；

4)建立红外光谱和仪器响应函数的数学模型：

$E(f,h)=\frac{1}{2}\left|\right|g-h\⊗f|{|}^2+\mathrm{\α}\Σ\mathrm{\ρ}(f^{\′})+\mathrm{\β}\Σ{\left|h_{\mathrm{\σ}}^{\′}\right|}^2$

上式中，第一项是数据项，其中E为响应函数，f为复原后的光谱，h为仪器响应函数，g表示退化的红外光谱；其余两项是正则化项，分别由正则化参数α和β来控制，其中ρ(·)为势函数， $\mathrm{\ρ}\left(t\right)=\left\{\begin{array}{cc}{t}^2& \left|t\right|\≤\mathrm{\μ}\\ 2\mathrm{\μ}\left|t\right|-{\mathrm{\μ}}^2& \left|t\right|>\mathrm{\μ}\end{array}\right.;$

其中μ是区分二次区域和线性区域的阈值参数；

5)采用交替最小值最优迭代算法求解红外光谱的数学模型，得到f和仪器响应函数h的宽度；；

6)判断迭代得到的光谱和仪器响应函数宽度是否收敛，若收敛，则输出红外光谱，否则返回步骤5)继续计算。

按上述方案，所述步骤3)中绝对中值估计器表达式δ，如下：

$\mathrm{\δ}=\frac{1.4826}{\sqrt{2}}median\left\{\left|g_i-g_{i-1}\right|,i=2,\mathrm{...},n\right\}$

其中，g表示退化的红外光谱。

按上述方案，所述正则化参数值初始化设置为α₀＝20δ，β₀∈[250,350]。

按上述方案，所述正则化参数值在每次迭代过程中按以下方式进行调节：

α＝α/λ₁，β＝β/λ₂；其中λ₁和λ₂为常数。

按上述方案，λ₁和λ₂的取值如下：λ₁＝1.01，λ₂＝1.02。

按上述方案，所述步骤6)中，收敛判断如下：模型求解完成后，判断f是否满足|σ_n+1-σ_n|＞d₁，||f_n+1-f_n||/||f_n||＞d₂；

其中d₁，d₂为非常小的正常数，d₁和d₂取值范围10^-9～10^-7；若满足，则判断为收敛。

本发明产生的有益效果是：

第一，提出了一种基于Huber-Markov约束的光谱反卷积的方法。对于红外光谱数据来讲，光谱噪声对原始曲线的谱带强度的影响是不相同的，在谱带强度变化平缓的波段受到噪声的影响较大，而在谱带强度变化剧烈的波段，受到噪声的影响较小。Huber函数根据谱带强度变化不同而采取不同强度的约束，在强度变化剧烈的地方采取较弱的约束，保存好光谱的结构；在强度变化平缓的地方采取较强的约束，抑制好噪声。这样就可以达到保存光谱结构的同时抑制光谱噪声。

第二，在光谱反卷积的约束中应用了Huber-Markov的自适应光滑性先验，能很好的保持光谱本身的结构，这是Huber-Markov先验在红外光谱数据反卷积中的首次得到应用。

第三，使用交替最小值最优的迭代算法来求解反卷积模型，该迭代算法克服了常规梯度降方法的不收敛的问题，具有很强的实用价值。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例的方法流程图；

图2是基于Huber-Markov先验的盲目光谱反卷积算法流图；

图3为红外光谱数据退化模型，其中：(a)退化的红外光谱g(b)原始的清晰的红外光谱f(c)仪器响应函数h(d)随机光谱噪声n；

图4为退化红外光谱分析图；其中：(a)原始的红外光谱，粗线代表陡峭区域，细线代表平坦区域(b)退化光谱，含有重叠谱带和噪声；

图5为不同的阈值μ情况下的Huber函数形态，可区分平坦区域和陡峭区域；

图6为原始的退化红外光谱曲线葡萄糖分子的(D(+)-Glucopy-ranose)的红外谱线(1500～10cm^-1)；

图7为归一化处理后的谱线[0,1]，计算光谱的一阶导数，并确定正则化参数的数值；

图8为红外光谱超分辨率复原实验结果，其中(a)为葡萄糖分子的(D(+)-Glucopyranose)的红外谱线(1500～10cm^-1)；(b)为使用FSD方法反卷积的结果，核的宽度是9.5cm^-1(c)HMSBD方法复原的结果；

图9为实测红外光谱实验的结果，其中(a)为波段从4000～1000cm^-1的化学物质脂肪族羧酸(Aliphatic,carboxylicacid,ether)的红外光谱(b)为FSD方法的反卷积结果，(c)为HMSBD方法的结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

通常，对于相干光和非相干光来说，光谱数据都是通过分光光度计来测量获得的。退化光谱可以被数学建模成实际光谱和仪器响应函数的卷积，再叠加上随机噪声。测量光谱和实际光谱之间的关系可以被表达为：

$g\left(v\right)=f\left(v\right)\⊗h\left(v\right)+n\left(v\right),$

其中g(v)是实际测量的光谱，而表示卷积操作；f(v)表示实际光谱，h(v)代表点扩散函数(PSF，也称为模糊核)，用来表达固有的现行变换函数和仪器扩散。通常来说h(v)的长度远小于f(v)。红外光谱退化示意图如图3所示，其光谱分析如图4所示，图中(a)为原始的红外光谱，粗线代表陡峭区域，细线代表平坦区域，(b)为退化光谱，含有重叠谱带和噪声。

对于光谱中使用Gauss-Markov先验(和Laplacian先验)，通常认为光谱中的高频成分会被去除掉。本发明中提出能够保存细节的Huber-Markov先验。这一先验能够有效地保存光谱中的细节。Huber-Markov先验和Gauss-Markov先验的唯一区别是势函数ρ(·)的不同，Huber函数定义如下：

$\mathrm{\ρ}\left(t\right)=\left\{\begin{array}{cc}{t}^2& \left|t\right|\≤\mathrm{\μ}\\ 2\mathrm{\μ}\left|t\right|-{\mathrm{\μ}}^2& \left|t\right|>\mathrm{\μ}\end{array}\right.$

其中μ是区分二次区域和线性区域的阈值参数，如图5所示。可以看出Gauss-Markov先验可以被当成Huber-Markov先验中μ趋向于+∞的特殊情况。

如图6显示了葡萄糖化合物的实测红外光谱曲线，曲线横坐标为波数，纵坐标为红外光谱强度，光谱分辨率为1cm^-1。谱线在400-600cm^-1波数范围内存在严重的谱带重叠现象，并且受到噪声的影响，该光谱可用来测试算法对光谱曲线的分辨率复原能力和对噪声的抑制能力。

参考图1，本发明方法具体流程为：

(1)将红外光谱的数据点记为求出光谱强度最大和最小的点的值，记为f_max和f_min，为了使算法适应不同类型的光谱曲线，对光谱作归一化处理，即

$f_i=\frac{{\hat{f}}_i-f_{\min }}{f_{\max }-f_{\min }}{f}_i\∈(0,1),i\∈(0,1,2,\mathrm{...}n)$

如图6所示，葡萄糖化合物的红外光谱在1500至400cm^-1的波数范围内，波长强度的最大值为34.而对于不同波数范围、不同物质其最大强度变化范围会有很大的不同，因此需对其进行归一化处理，对图6进行归一化处理后的结果如图7所示。

(2)求归一化后光谱强度f_i的一阶导数，f_i′＝(f_i+1-f_i)/2，相邻的点做差分计算；

(3)根据上式光谱一阶导数的强度值，使用绝对中值估计器估计光谱的噪声水平(步骤6所示的公式)，确定正则化参数的数值；计算出δ＝0.0015，设置正则化参数为α₀＝0.03，β₀＝300。初始化仪器响应函数h＝G_σ，其中σ＝1；

(4)根据前面对红外光谱特性的分析，本发明提出的数学模型：

$E(f,h)=\frac{1}{2}\left|\right|g-h\⊗f|{|}^2+\mathrm{\α}\Σ\mathrm{\ρ}(f^{\′})+\mathrm{\β}\Σ|h_{\mathrm{\σ}}^{\′}{|}^2$

上式中，第一项是数据项，其余两项是正则化项，分别由正则化参数α和β来控制。下面使用交替最小值最优迭代方法求解该数学模型，解得到红外光谱分辨率；

(5)固定h而优化f。丢弃掉常数项之后数学模型E(f,h)可以简写成E(f)，如下式：

$E\left(f\right)=\frac{1}{2}\left|\right|g-h\⊗f|{|}^2+\mathrm{\α}\Σ\mathrm{\ρ}\left(f^{\′}\right)$

该式在实数域IR上是凸的，而且如果解存在，解将是唯一的并且可以用梯度下降法来获得。该等式的欧拉-拉格朗日(E-L)方程可表达为：

$\frac{\mathrm{\δ}E}{\mathrm{\δ}f}=\underset{i\∈L}{\Σ}\left(\frac{\∂F}{\∂f}-\frac{d}{dx}\left(\frac{\∂F}{\∂f^{\′}}\right)\right).$

上式中对应的两项的梯度分别表示为：

$\frac{\∂F}{\∂f}=(f\⊗h_{\mathrm{\σ}}-g)\⊗h_{\mathrm{\σ}}(-v)$

和

$\frac{d}{dx}\left(\frac{\∂F}{\∂f^{\′}}\right)=\left\{\begin{array}{cc}2f^{\′\′}& \left|f^{\′}\right|\≤\mathrm{\μ}\\ 0& \left|f^{\′}\right|>\mathrm{\μ}\end{array}\right..$

因此，欧拉方程可以重写成：

$\frac{\mathrm{\δ}E}{\mathrm{\δ}f}=\left\{\begin{array}{cc}\left(f\⊗h_{\mathrm{\σ}}-g\right)\⊗h_{\mathrm{\σ}}\left(-v\right)-2{\mathrm{\αf}}^{\′\′}& \left|f^{\′}\right|\≤\mathrm{\μ}\\ \left(f\⊗h_{\mathrm{\σ}}-g\right)\⊗h_{\mathrm{\σ}}\left(-v\right)& \left|f^{\′}\right|>\mathrm{\μ}.\end{array}\right.$

连续逼近迭代表达式如下：

${\hat{f}}_{n+1}={\hat{f}}_n+t_n\left(-\frac{\mathrm{\δ}E}{{\mathrm{\δf}}_n}\right)$

n表示迭代次数，t_n代表时间步长。如果t_n取得太小，迭代速度将会十分缓慢；反之算法不收敛。使用二阶Taylor级数逼近目标函数，二次步长的逼近表达式为：

$t_n=\frac{{(\▿E(f_n\left)\right)}^T\▿E\left(f_n\right)}{(\▿E(f_n\left)\right)\left({\▿}^{2}E\right(f_n\left)\right)\▿E\left(f_n\right)}$

其中和分别是目标函数E(f)的一阶矩阵和Hessian矩阵。

(6)固定f来计算最优的仪器响应函数的宽度σ，数学模型可简化为

$E\left(h\right)=\frac{1}{2}\left|\right|h_{\mathrm{\σ}}\⊗f-g|{|}^2+\mathrm{\β}\Σ|h_{\mathrm{\σ}}^{\′}{|}^2$

相对于标量参数σ，对其求偏导数

$\frac{\∂E}{\∂\mathrm{\σ}}=\underset{i\∈L}{\Σ}(h_{\mathrm{\σ}}\⊗f-g){\left(\frac{\∂h_{\mathrm{\σ}}}{\∂\mathrm{\σ}}\⊗f\right)}_i+\mathrm{\β}\underset{i\∈L}{\Σ}{\left(\frac{\∂}{\∂\mathrm{\σ}}|{h_{\mathrm{\σ}}}^{\′}{|}^2\right)}_i=0,$

上式中

$\frac{\∂h_{\mathrm{\σ}}}{\∂\mathrm{\σ}}=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}^2}\exp (-\frac{v^2}{2{\mathrm{\σ}}^2})\·\left(\frac{v^2}{{\mathrm{\σ}}^2}-1\right)$

$h_{\mathrm{\σ}}^{\′}=\frac{\∂h_{\mathrm{\σ}}}{\∂v}=-\frac{v}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}^3}\exp (-\frac{v^2}{2{\mathrm{\σ}}^2}),$

$\frac{\∂}{\∂\mathrm{\σ}}|h_{\mathrm{\σ}}^{\′}{|}^2=\frac{v^2}{{\mathrm{\π\σ}}^7}\exp \left(-\frac{v^2}{{\mathrm{\σ}}^2}\right)\left(\frac{v^2}{{\mathrm{\σ}}^2}-3\right).$

上式中σ通过二分法直接求解。

(7)f，h求解完成后，更新h，判断f是否满足|σ_n+1-σ_n|＞d₁,||f_n+1-f_n||/||f_n||＞d₂，其中d₁，d₂为非常小的正常数，取值范围为[10^-9，10^-7]，若不满足条件，则返回步骤(5)继续求解，若满足条件，则输出高分辨率红外光谱f^k+1；其流程图如图2所示。

本发明中，两个正则化参数(α和β)是可调节的。α和β对应于数学模型中光谱与核先验的概率参数，而且它们的值可以在优化迭代过程中被改变。我们先设置了α＝α/λ₁和β＝β/λ₂，然后经过每一次迭代，α和β的值分别被除以λ₁和λ₂。我们通常设置λ₁＝1.01和λ₂＝1.02，用来在迭代后期减少对光谱结构的影响而增加与原始光谱的相似性。

为了验证本发明算法的有效性，本发明做了对比实验，对比方法有傅里叶自反卷积(FSD)(见文献：J.K.Kauppinen,D.J.Moffatt,H.H.Mantsch,andD.G.Cameron,“FourierSelf-Deconvolution:AMethodforResolvingIntrinsicallyOverlappedBands,”Appl.Spectrosc.,vol.35,no.3,pp.271-276,1981.)、半盲目反卷积方法(SBD)(见文献：L.Yan,H.Liu,S.Zhong,andH.Fang,“Semi-BlindSpectralDeconvolutionwithAdaptiveTikhonovRegularization,”AppliedSpectroscopy,vol.66,no.11,pp.1334-1346,2012.)。本发明中，简称本发明的方法名称为HMSBD(Huber-MarkovSpectralBlindDeconvolution)，实验结果曲线如图8和图9所示。本发明的方法在细节保存和噪声抑制两方面都优于其它方法。在提高光谱分辨率方面，我们计算了光谱的半高宽比，显示在表1中。根据半高宽比提高的程度，本发明的方法显著优于其它方法，光谱分辨率平均提高了约30％。

表1三种算法的半高宽比

待处理红外光谱	原始	FSD	SBD	MOBSD
					红外光谱1	1.75	1.89	2.33	3.11
红外光谱2	1.40	1.46	2.15	2.87
					红外光谱3	1.35	1.48	1.69	2.32

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (6)

1.一种盲目反卷积红外光谱图超分辨率复原方法，其特征在于，包括以下步骤：

1)采集待处理的红外光谱，对待处理的红外光谱数据作归一化处理，其中i表示波数点，n表示波数的点数；

2)对归一化后光谱强度f_i进行一阶求导，f′_i＝(f_i+1-f_i)/2；

3)根据光谱一阶导数强度值f′_i，使用绝对中值估计器初步计算光谱噪声水平，并由此确定正则化参数的数值；

4)建立红外光谱和仪器响应函数的数学模型：

$E(f,h)=\frac{1}{2}\left|\right|g-h\⊗f|{|}^2+\mathrm{\α}\mathrm{\Σ}\mathrm{\ρ}\left(f^{\′}\right)+\mathrm{\β}\mathrm{\Σ}|h_{\mathrm{\σ}}^{\′}{|}^2$

上式中，第一项是数据项，其中f为复原后的光谱，h为仪器响应函数，g表示退化的红外光谱；其余两项是正则化项，分别由正则化参数α和β来控制，

其中ρ(·)为势函数， $\mathrm{\ρ}\left(t\right)=\left\{\begin{array}{cc}{t}^2& \left|t\right|\≤\mathrm{\μ}\\ 2\mathrm{\μ}\left|t\right|-{\mathrm{\μ}}^2& \left|t\right|>\mathrm{\μ}\end{array}\right.;$

其中μ是区分二次区域和线性区域的阈值参数；

5)采用交替最小值最优迭代算法求解红外光谱的数学模型，得到f和仪器响应函数h的宽度；

6)判断迭代得到的光谱和仪器响应函数宽度是否收敛，若收敛，则输出红外光谱，否则返回步骤5)继续计算。

2.根据权利要求1所述的复原方法，其特征在于，所述步骤3)中绝对中值估计器表达式δ，如下：

$\mathrm{\δ}=\frac{1.4826}{\sqrt{2}}median\left\{|g_i-g_{i-1}|,i=2,...,n\right\}$

其中，g表示退化的红外光谱。

3.根据权利要求1所述的复原方法，其特征在于，所述正则化参数值初始化设置为α₀＝20δ，β₀∈[250,350]。

4.根据权利要求1所述的复原方法，其特征在于，所述正则化参数值在每次迭代过程中按以下方式进行调节：

α＝α/λ₁，β＝β/λ₂；其中λ₁和λ₂为常数。

5.根据权利要求4所述的复原方法，其特征在于，λ₁和λ₂的取值如下：λ₁＝1.01，λ₂＝1.02。

6.根据权利要求1所述的复原方法，其特征在于，所述步骤6)中，收敛判断如下：模型求解完成后，判断f是否满足|σ_n+1-σ_n|＞d₁，||f_n+1-f_n||/||f_n||＞d₂；其中d₁，d₂为非常小的正常数，d₁和d₂取值范围10^-9～10^-7；若满足，则判断为收敛。