一种基于自适应正则化的光谱分辨率增强方法

本发明涉及一种基于自适应正则化的光谱分辨率增强方法，属于光谱技术 与信号处理交叉的领域。

光谱分析作为一种快速、无损、多组分同时分析的技术，在许多领域中被 广泛地应用，如拉曼光谱广泛地用于食品质量监测、化学品快速检测、材料分 析和生物医学等领域。然而，光谱数据在光谱仪测量过程中，经常会受到谱线 的自然展宽、碰撞展宽、探测器和电路低通特性等因素的影响，会导致光谱数 据强噪声和子带重叠现象。这些干扰使得光谱的分辨率下降，给光谱的分析和 应用带来了困难，极大的限制了光谱的应用。为最大程度地反映物质本来的光 谱特征，尽可能的通过一定的光谱信号处理方法和手段，使得光谱的分辨率得 到增强。

目前，对退化光谱进行分辨率增强的方法有很多，主要可分为两大类：曲 线拟合和光谱反卷积。曲线拟合的方法认为光谱曲线是通过高斯函数和洛伦兹 函数的组合而成，即实测光谱能通过这些函数来拟合得到。这种方法能够很好 地解决子带重叠的问题，即能拟合出重叠的波峰，但是对噪声极其敏感。光谱 反卷积的方法主要是通过不同的约束函数来进行光谱分辨率增强，这些方法有 高阶统计、Tikhonov正则化、同态滤波等。在这些方法中，当函数的约束能力 较强时，能够很好地抑制光谱噪声，但过强的约束会导致体现光谱特征的细节 信息(如谱峰)也会被平滑，从而不能很好地增强光谱分辨率。当函数的约束 能力较弱时，虽然能够将退化的细节较好的恢复出来，但抑制噪声能力会减弱。 因此，这类方法在光谱的分辨率增强和噪声抑制能力上很难达到平衡。而本发 明能够很好地解决上面的问题。

本发明目的在于提供了一种基于自适应正则化的反卷积方法，该方法能够 使光谱分辨率增强，解决了由光谱仪器测得的光谱数据经常会受到随机噪声和 仪器响应函数误差的影响而导致光谱退化、分辨率降低的问题。

本发明解决其技术问题所采取的技术方案是：一种基于自适应正则化的光 谱分辨率增强方法，该方法基于自适应正则化的反卷积模型，用于光谱分辨率 增强，能够为光谱在质量检测和材料分析等方面的应用提供技术支撑。

方法流程：

步骤1：对离散化的光谱数据 $\hat{f}=(\widehat{f_1},\widehat{f_2},\widehat{f_3}...\widehat{f_{i-1}},\widehat{f_i},\widehat{f_{i+1}}...\widehat{f_{N-2}},\widehat{f_{N-1}},\widehat{f_N})$ 进行归 一化预处理：

$f_i=\frac{\widehat{f_i}-f_{\min }}{f_{\max }-f_{\min }},f_i\in \left(\mathrm{0,1}\right),i\in (\mathrm{0,1},...,N)---\left(1\right)$

其中，fmax和fmin分别为光谱强度最大和最小值，得到的归一化后的光谱数据为 f＝(f1,f2,f3...fi-1,fi,fi+1...fN-2,fN-1,fN)。

步骤2：利用中心差分法求归一化后光谱强度fi的一阶导数(假设分辨率为 1)：

${f_i}^{\prime }=\frac{f_{i+1}-f_{i-1}}{2},i\in (\mathrm{0,1},...,N)---\left(2\right)$

步骤3：设定仪器响应函数hσ(x)，本发明选择带参数的高斯线形函数作为 模糊核函数：

$h_{\sigma }\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{x^2}{2{\sigma }^2})---\left(3\right)$

这里σ为高斯核参数，控制卷积核的宽度(高斯卷积核的半高宽的一半宽度为
则

$\left|h_{\sigma }^{\prime }\left(v\right)\right|=\frac{\left|v\right|}{\sqrt{2\pi }{\sigma }^3}\exp (-\frac{v^2}{2{\sigma }^2})---\left(4\right)$

步骤4：本发明提出的自适应光谱反卷积的模型构建为：

$E(f,h_{\sigma })=\frac{1}{2}{\left|\right|f*h_{\sigma }-g\left|\right|}^2+\mathrm{\alpha \Sigma }\ln (1+w{\left|f^{\prime }\right|}^2)+\mathrm{\beta \Sigma }\left|h_{\sigma }^{\prime }\right|---\left(5\right)$

其中，“||●||”表示求L2范数，α和β分别为控制谱的平滑性和模糊核函数的核 宽度参数σ的正则化参数。w为用于保留光谱细节的自适应系数

$w_i=\frac{k^2}{k^2+{f_i}^2},i\in (\mathrm{0,1},...,N)---\left(6\right)$

其中，k为一常数。可以看出：在光谱的平坦区，f较小，w相对较大，使 得正则项加于该区的强度较强，可以很好的平滑噪声；在光谱的峰值附近区域， f相对较大，w相对较小，使得正则项加于该区的强度较弱，从而很好的保存 谱峰处的细节信息。

步骤5：本发明将求解复原光谱及核参数的问题转化为求解最小化问题：

$(\hat{f},\widehat{h_{\sigma }})=\arg \underset{f,h_{\sigma }}{\min }E(f,h_{\sigma })---\left(7\right)$

并采用交替最小化方法进行求解，即

$\left\{\begin{array}{cc}\hat{f}=\arg \underset{f}{\min }{E}_1\left(f\right)& \mathrm{fix}{h}_{\sigma }\\ \widehat{h_{\sigma }}=\arg \underset{h_{\sigma }}{\min }{E}_2\left(h_{\sigma }\right)& \mathrm{fix\; f}\end{array}\right.---\left(8\right)$

其中，

$E_1\left(f\right)=\frac{1}{2}{\left|\right|f*h_{\sigma }-g\left|\right|}^2+\mathrm{\alpha \Sigma }\ln (1+w{\left|f^{\prime }\right|}^2)---\left(9\right)$

$E_2\left(f\right)=\frac{1}{2}{\left|\right|f*h_{\sigma }-g\left|\right|}^2+\mathrm{\beta \Sigma }\left|h_{\sigma }^{\prime }\right|---\left(10\right)$

具体步骤如下：

①设置参数α,β,k的初始值，给定最大迭代次数nmax；

②初始化n＝0，f0＝g，σ＝σ0＝1，σk≥1；

③迭代求解：

(i)令σk＝σ，n＝n+1；

(ii)令hσ不变，最小化E1(f)求解fn+1：

利用泛函求解式(9)的最小化问题。根据变分原理，有

$\delta E_1=\Sigma \left[(\frac{\partial F}{\partial f}-\frac{d}{\mathrm{dv}}\left(\frac{\partial F}{\partial f^{\prime }}\right))\mathrm{\delta f}\right]=0---\left(11\right)$

其中，

$\frac{\partial F}{\partial f}=(f*h_{\sigma }-g)*h_{\sigma }(-v)-2\alpha \frac{k^{2}f{\left|f^{\prime }\right|}^2}{(k^2+{\left|f\right|}^2+k^2{\left|f^{\prime }\right|}^2)(k^2+{\left|f\right|}^2)}---\left(12\right)$

$\frac{d}{\mathrm{dv}}\left(\frac{\partial F}{\partial f^{\prime }}\right)=2\alpha k^2\frac{k^2{f}^{\prime \prime }+{\left|f\right|}^2{f}^{\prime \prime }-2f{\left|f^{\prime }\right|}^2-k^2{\left|f^{\prime }\right|}^2{f}^{\prime \prime }}{{(k^2+{\left|f\right|}^2+k^2{\left|f^{\prime }\right|}^2)}^2}---\left(13\right)$

f″＝fi+1-2fi+fi-1    (14)

本发明利用梯度下降法对式(11)进行求解，即

$f^{n+1}=f^n+\mathrm{\Delta t}(-\frac{\delta E_1}{\mathrm{\delta f}})---\left(15\right)$

其中，n为迭代次数，Δt通常设为0.5～0.8，为步长因子。

(iii)令f＝fⁿ⁺¹不变，由于模糊核函数h_σ采用参数化形式，其主要取决于
参数σ，因此最小化E₂(h_σ)可通过求解得到，即

$\frac{\partial E_2}{\partial \sigma }=\Sigma [(f*h_{\sigma }-g)(\frac{\partial h_{\sigma }}{\partial \sigma }*f)+\beta \frac{\partial \left|h_{\sigma }^{\prime }\right|}{\partial \sigma }]=0---\left(16\right)$

本发明采用二分法求解式(16)得到高斯核参数σ：设定一个二分区间σmin 和σmax，给定迭代次数Nσ，并假设a＝σmin，b＝σmax；令 $I=\Sigma [(f*h_{\sigma }-g)(\frac{\partial h_{\sigma }}{\partial \sigma }*f)+\beta \frac{\partial \left|h_{\sigma }^{\prime }\right|}{\partial \sigma }];$

当nσ＜Nσ时，进行如下迭代：

${\sigma }_{\min }=\frac{a+b}{2},n_{\sigma }=n_{\sigma }+1;$

判断|b-σ_mid|≤ε，若是，则停止迭代；否则，将b和σ_mid带入求得对应的I_b
和；判断，若是，则a＝σ_mid；否则，b＝σ_mid；

(iv)判断收敛条件：如果|σk-σ|＞ε(ε为设定的一个小量)且n＜nmax，则： 重复(i)～(iii)；否则，停止迭代，输出分辨率增强后的光谱和模糊核函数的 估计结果。

本发明用于光谱分辨率增强，其既能抑制噪声的同时增强光谱分辨率，又 能自适应估计模糊核函数宽度。其技术效果体现在：第一，充分利用了噪声对 光谱强度受的影响远小于对光谱强度一阶导数的影响的特点，对光谱信息的正 则项施加了自适应权系数，更好的保存光谱峰值；第二，利用光谱的先验知识， 将模糊核函数构造成参数化形式，缓解了盲卷积计算量大和非盲反卷积鲁棒性 弱的缺点；第三，利用交替迭代最小化方法迭代求解，并分别用二分法求解仪 器响应函数和用梯度下降法求解分辨率增强后的光谱曲线，使得算法具有很好 的鲁棒性。本发明对光谱在质量检测和材料分析等方面提供了技术支撑，具有 较强的实用价值。

图1为本发明的方法流程图。

图2为真实光谱归一化曲线图。

图3为模糊核函数曲线图。

图4为模糊化后的无噪声光谱曲线图。

图5为模糊化后带噪声的光谱曲线图。

图6为无噪声模糊光谱恢复结果曲线图。

图7为带噪声模糊光谱恢复结果曲线图。

图8为对无噪声模糊光谱恢复的模糊核自适应估计结果曲线图。

图9为对待噪声模糊光谱恢复的模糊核自适应估计结果曲线图。

图10为L(+)-Arabinofuranose(C5H10O5)的实测拉曼谱曲线(从1400cm-1到 900cm-1)曲线图。

图11为L(+)-Arabinofuranose(C5H10O5)实测拉曼谱的分辨率增强结果曲线 图。

下面结合附图和实例对本发明作进一步详细的说明。

如图1所示，本发明提供了一种基于自适应正则化的光谱分辨率增强方法， 该方法是通过利用光谱的先验知识构造自适应正则项和参数化模糊核函数，建 立反卷积模型，包括如下步骤：

步骤1：对离散化的光谱数据 $\hat{f}=(\widehat{f_1},\widehat{f_2},\widehat{f_3}...\widehat{f_{i-1}},\widehat{f_i},\widehat{f_{i+1}}...\widehat{f_{N-2}},\widehat{f_{N-1}},\widehat{f_N})$ 进行归 一化预处理得到的归一化后的光谱数据f＝(f1,f2,f3...fi-1,fi,fi+1...fN-2,fN-1,fN)， 并利用中心差分法求归一化后光谱强度fi的一阶导数fi′；

步骤2：为了平滑噪声，构造关于光谱数据的正则项

αΣln(1+|f′|2)    (17)

式(17)所示的正则项对f′较大的区域施加较大的正则化力，而对f′较小 的区域施加较小的正则化力；而光谱细节处的f′较大，直接用式(17)所示的 正则项在平滑噪声的同时，也会平滑光谱细节；

步骤3：为了进一步保留光谱细节，提出用于保留光谱细节的自适应系数

$w_i=\frac{k^2}{k^2+{f_i}^2},i\in (\mathrm{0,1},...,N)---\left(18\right)$

其中，k为一常数；

步骤4：设定仪器响应函数hσ(x)，本发明选择带参数的高斯线形函数作为模 糊核函数：

$h_{\sigma }\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }\exp (-\frac{x^2}{2{\sigma }^2})---\left(19\right)$

其中，σ为高斯核参数，控制卷积核的宽度，高斯卷积核的半高宽的一半宽
度为

步骤5：本发明提出的基于自适应正则化的光谱分辨率增强模型构建为：

$E(f,h_{\sigma })=\frac{1}{2}{\left|\right|f*h_{\sigma }-g\left|\right|}^2+\mathrm{\alpha \Sigma }\ln (1+w{\left|f^{\prime }\right|}^2)+\mathrm{\beta \Sigma }\left|h_{\sigma }^{\prime }\right|---\left(20\right)$

其中，“||●||”表示求L2范数，α和β分别为控制谱的平滑性和模糊核函数的 核宽度参数σ的正则化参数；

步骤6：本发明将求解复原光谱及核参数的问题转化为求解最小化问题：

$(\hat{f},\widehat{h_{\sigma }})=\arg \underset{f,h_{\sigma }}{\min }E(f,h_{\sigma })---\left(21\right)$

并采用交替最小化方法进行求解，即

$\left\{\begin{array}{cc}\hat{f}=\arg \underset{f}{\min }{E}_1\left(f\right)& \mathrm{fix}{h}_{\sigma }\\ \widehat{h_{\sigma }}=\arg \underset{h_{\sigma }}{\min }{E}_2\left(h_{\sigma }\right)& \mathrm{fix\; f}\end{array}\right.---\left(22\right)$

其中，

$E_1\left(f\right)=\frac{1}{2}{\left|\right|f*h_{\sigma }-g\left|\right|}^2+\mathrm{\alpha \Sigma }\ln (1+w{\left|f^{\prime }\right|}^2)---\left(23\right)$

$E_2\left(f\right)=\frac{1}{2}{\left|\right|f*h_{\sigma }-g\left|\right|}^2+\mathrm{\beta \Sigma }\left|h_{\sigma }^{\prime }\right|---\left(24\right)$

本发明所述方法对上述式(23)和(24)的求解过程包括：

(A)令hσ不变，利用泛函最小化E1(f)求解fn+1。根据变分原理，得到：

$\delta E_1=\Sigma \left[(\frac{\partial F}{\partial f}-\frac{d}{\mathrm{dv}}\left(\frac{\partial F}{\partial f^{\prime }}\right))\mathrm{\delta f}\right]=0---\left(25\right)$

其中，

$\frac{\partial F}{\partial f}=(f*h_{\sigma }-g)*h_{\sigma }(-v)-2\alpha \frac{k^{2}f{\left|f^{\prime }\right|}^2}{(k^2+{\left|f\right|}^2+k^2{\left|f^{\prime }\right|}^2)(k^2+{\left|f\right|}^2)}---\left(26\right)$

$\frac{d}{\mathrm{dv}}\left(\frac{\partial F}{\partial f^{\prime }}\right)=2\alpha k^2\frac{k^2{f}^{\prime \prime }+{\left|f\right|}^2{f}^{\prime \prime }-2f{\left|f^{\prime }\right|}^2-k^2{\left|f^{\prime }\right|}^2{f}^{\prime \prime }}{{(k^2+{\left|f\right|}^2+k^2{\left|f^{\prime }\right|}^2)}^2}---\left(27\right)$

f″＝fi+1-2fi+fi-1    (28)

本发明所述方法利用梯度下降法对式(25)进行求解，即

$f^{n+1}=f^n+\mathrm{\Delta t}(-\frac{\delta E_1}{\mathrm{\delta f}})---\left(29\right)$

其中，n为迭代次数，Δt通常设为0.5～0.8，为步长因子；

(B)令f＝fⁿ⁺¹不变，由于模糊核函数h_σ采用参数化形式，其取决于参数σ，
因此最小化E₂(hσ)可通过求解得到，即

$\frac{\partial E_2}{\partial \sigma }=\Sigma [(f*h_{\sigma }-g)(\frac{\partial h_{\sigma }}{\partial \sigma }*f)+\beta \frac{\partial \left|h_{\sigma }^{\prime }\right|}{\partial \sigma }]=0---\left(30\right)$

本发明所述方法是采用二分法求解式(30)得到高斯核σ。

通常观测得到的光谱可描述成真实光谱和测量仪器响应函数(又称为模糊 核函数)的卷积，附加噪声组成，用数学表达式可表示成：

g(v)＝f(v)*h(v)+n(v)    (31)

其中，g(v)为观测到的光谱数据，f(v)为真实光谱数据，h(v)为仪器响应函数， “*”表示卷积算子，即f(v)*h(v)＝Σkfkh(v-k)，n(v)是测量时引入的噪声， 通常为高斯白噪声，v为波数(wavenumber)，单位为cm-1。为了验证本发明对 光谱分辨率增强的有效性，首先将已知的某红外光谱进行模糊化，然后用本发 明提供的方法进行恢复，最后将恢复结果与真实光谱进行比较。

图2给出了某真实红外光谱归一化后的曲线图，其分辨率为1cm-1，将其与 图3所示的模糊核进行卷积，得到如图4所示的模糊化结果；在将其加上高斯 白噪声，得到如图5所示的带噪声的模糊化结果。比较图2和图4、图5可以发 现，模糊化后的光谱图相比于真实光谱分辨率降低了，特别是2924cm-1和2852 cm-1处的两个谱峰重叠在一起，3364cm-1到2284cm-1波数范围内，以及1464cm-1 到1338cm-1波数范围内也都存在严重的变形和重叠。

将图4和图5作为待处理的观测光谱数据，按照以下步骤进行处理：

(1)读入光谱数据，如果没有数据输入，则转至步骤(6)；

(2)初始化参数：α＝0.2，β＝0.65，k＝0.1，给定最大迭代次数nmax＝500； 并令初始迭代次数n＝0，复原光谱的初始值f0＝g，模糊核参数的初始值为 σ＝σ0＝1，σk≥1；

(3)σ＝σn，σk＝σ利用梯度下降法求解fn+1；

(4)令f＝fn+1不变，采用二分法求解高斯核参数σn+1；

(5)n＝n+1，判断收敛条件：如果|σk-σ|＞ε且n≤nmax，则重复上述步骤 (3)和(4)；否则，停止迭代，输出复原光谱和模糊核函数的估计结果。

(6)算法结束。

利用上述处理步骤，对应于图4和图5，得到复原光谱曲线如图6和图7所 示，自适应估计模糊核函数的结果如图8和图9所示。将图6和图7与图2进 行比较、图8和图9与图3进行比较可以看出，用本发明提供的方法能够对模 糊化后的光谱进行很好的复原，且能准确的估计模糊核函数。

为了进一步量化说明本发明提供方法的复原效果，本发明分别将图6、图7
与图2进行相似性比较，求得表征相关相似度的相关系数
<math> <mrow> <mrow> <mo>(</mo> <mi>CC</mi> <mo>=</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>L</mi> </msubsup> <mrow> <mo>(</mo> <mover> <msub> <mi>f</mi> <mi>i</mi> </msub> <mo>^</mo> </mover> <mo>-</mo> <mover> <mover> <mi>f</mi> <mo>^</mo> </mover> <mo>&OverBar;</mo> </mover> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>f</mi> <mo>&OverBar;</mo> </mover> <mo>)</mo> </mrow> <mo>/</mo> <msqrt> <mo>[</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </msubsup> <msup> <mrow> <mo>(</mo> <mover> <msub> <mi>f</mi> <mi>i</mi> </msub> <mo>^</mo> </mover> <mo>-</mo> <mover> <mover> <mi>f</mi> <mo>^</mo> </mover> <mo>&OverBar;</mo> </mover> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>]</mo> <mo>[</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </msubsup> <msup> <mrow> <mo>(</mo> <msub> <mi>f</mi> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>f</mi> <mo>&OverBar;</mo> </mover> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>]</mo> </msqrt> </mrow> <mo>,</mo> </mrow> </math> 其中，
分别为复原光谱与原始光谱的均值)分别为0.9969和
0.9953。可以看出，相关系数基本都趋近于1，说明复原的光谱与真实光谱是非
常相近的。

图10给出了从物质L(+)-Arabinofuranose(C5H10O5)的实测拉曼谱曲线截取 的从1400cm-1到900cm-1的结构，图11给出了用本发明提出的基于自适应正 则化的光谱分辨率增强方法对其图10所示光谱曲线的增强结果。可以明显看出， 图11所示的光谱的分辨率明显高于图10所示的光谱分辨率。

以上所述仅为本发明的一个实施例而已，并不用以限制本发明，凡在本发明 的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的 保护范围之内。