CN105205788B - 一种针对高通量基因测序图像的去噪方法 - Google Patents

Abstract

本发明提出了一种基于à trous小波阈值去噪方法，应用于高通量基因测序图像去噪。本发明主要是对基于硬阈值的各向同性非抽取离散小波算法的改进，具体操作是使用特定的小波对高通量基因测序的图像进行小波分解，针对每一层小波系数，提出使用l₁范数计算全局阈值。通过全局阈值和每一层小波系数构建估计小波系数表达式，最后使用小波重构算法得到去噪后的图像。本发明的方法相对于现有技术中具有去噪性能更好的优点，并且通过几组实验数据对比，证明本发明的合理性和鲁棒性。

Description

一种针对高通量基因测序图像的去噪方法

技术领域

本发明属于数字图像降噪领域，具体涉及一种针对高通量基因测序图像的去噪方法。

背景技术

高通量基因测序图像承载着丰富的人类基因信息，极高的清晰度要求已经成为衡量实验成功的重要环节。但是由于高通量基因测序图像在获取或者传输的过程中，会生成多种类型的噪声，这些噪声影响后续的图像处理、点检测、以及Base Calling等操作，因此，对图像进行去噪具有重大意义。而主要噪声包括如下两类：①高通量基因测序的原始图像通过CCD相机获得，当光通过传感器产生信号电荷的过程视为随机过程，而在单位时间内，电荷数目在平均值上小的波动被认为是泊松噪声。②待测序的图像含有的噪声具有随机性或者概率性，在实际应用中，经常会被建模成为均值为零的白噪声，而图像中碱基灰度分布满足高斯分布，因此这种白噪声就会是高斯白噪声。所以，高通量基因测序的图像含有的噪声包括高斯噪声、泊松噪声。

通过CCD相机获得的高通量基因测序图像中，由于待测序列的每个碱基都受到荧光蛋白的标记，在图像中显示成为由几个像素组成的亮点。所以图像是由多个大小不同的亮斑组成，并且图像极具纹理密度多样性的特性。由于获得图片质量受到设备的限制，原始图片的信噪比很低。目前针对待测序图像的去噪算法分为两类：空间域的图像去噪算法和频率域的图像去噪算法。空间域的图像去噪算法包括：高斯平滑滤波器，顶帽变换滤波器等；频率域的图像去噪算法包括：基于硬阈值的小波阈值收缩法，基于软阈值的小波阈值收缩法等。

(1)高斯平滑滤波

若对原图I使用高斯核G_σ进行图像平滑去噪，滤波后的图像J应表示为：

公式中*代表卷积符号。针对图像中不相干的噪声，这种平滑去噪方式与滤波器的选取有关，经选择的滤波器可以使得去噪后的图像信噪比最大。这是因为点扩散函数(Point Spread Function，PSF)修改了每个亮点的像素分布，进而能够通过高斯平滑得到更好的去噪效果。

(2)顶帽变换滤波

Yoshitaka Kimori等人提出将改进的顶帽变换算法用于高通量基因测序图像的去噪和点识别。通过将原图按照逆时针旋转N次相同的角度得到N幅图像，每幅图像使用线性结构元素进行开运算，然后将经过处理后的图像按顺时针的方向还原成原图像，选取N幅图像相同位置最大的灰度值构成开运算后的图像，最后用原图减去经过开运算后的图像得到实验结果。该方法通过选取大小合适的线性结构元素来抑制噪声，并能有效的进行点检测。

(3)基于硬阈值的小波阈值收缩法

Donoho和Johnstone等人提出小波收缩法，小波变换主要是通过获得少量较大数值的小波系数，进而获得真实信号中较大的能量，并且将小波系数中因为噪声引起的较小数值丢弃。因此小波收缩函数具有两个特征：1.舍弃数值小的小波系数；2.保留数值大的小波系数。小波收缩法分为阈值去噪法和比例去噪法，而阈值去噪法是其中较常用的方法。基于硬阈值的小波阈值收缩法是将一副图像进行小波分解后，得到不同频率的小波系数，将得到的小波系数与设计的阈值进行对比，通过以下公式得到估计小波系数，对于大于阈值的小波系数保留，对于小于阈值的小波系数归零。最后进行图像重构，得到去噪图像。

(4)基于软阈值的小波阈值收缩法

基于软阈值的小波阈值收缩法是将一副图像进行小波分解后，得到不同频率的小波系数，将得到的小波系数与设计的阈值进行对比，通过以下公式得到估计小波系数，对于大于阈值的小波系数要减去阈值系数，对于小于阈值的小波系数归零。最后进行图像重构，得到去噪后的图像。

以下是几种常用的阈值计算公式：

①VisuShrink阈值

D.L.Dononho和L.M.Johnstone在1994年提出Visushrink阈值，Visushrink阈值为：

其中σ代表噪声标准差，N代表信号的长度。

②Gaussian分布的置信区间阈值

零均值的Gaussian分布变量大部分都会落在[-3σ,3σ]，落在这个区间外的概率很小。因此，通过选择λ＝3σ～4σ，小波系数的绝对值小于阈值会被认为是噪声，而大于阈值的小波系数被认为主要由信号系数组成。

由于噪声的多样性，高斯滤波器大部分只能去除一种类型的噪声，而且在消除噪声的同时，图像噪声的边缘会变模糊。而顶帽变换滤波器该方法通过选取大小合适的线性结构元素来抑制噪声，并能有效的进行点检测。但是由于待测序碱基的大小不是固定的，所以碱基的大小影响线性结构元素的选取，进而影响顶帽变换滤波器的性能。而小波硬阈值收缩法可以得到较好的局部特征，但是由于估计小波系数的分布含有两个断点，会造成振铃、伪吉布斯等视觉失真现象。而且该算法会随着数据的微小变换而变换，所以该算法会产生较大的方差和不稳定性现象。而小波软阈值收缩算法的估计系数虽然整体的连续性好，使得去噪效果相对平滑。但是该算法仍然存在缺点，会随着大的小波阈值的收缩，产生较大的偏差。并且估计系数与真实系数之间存在恒定的误差，造成重构图像出现不必要的误差。

发明内容

本发明的目的在于提供一种用于去除高通量测序图像噪声的方法，旨在解决现有技术对高通量基因测序图像进行去噪操作时存在的问题，比如去噪效果不好，去噪后的图片存在失真现象。为了实现目的，本发明提供了一种基于新的小波阈值收缩的图像去噪方法，该方法是通过在频率域内实现对图像进行去噪处理，相对于现有技术中具有去噪性能更好的优点。

本发明具体通过如下技术方案实现：

一种针对高通量基因测序图像的去噪方法，其包括以下步骤：

(1)对测序图像使用小波函数进行小波分解，得到每一层小波系数ω_i；

(2)针对每一层小波系数，计算出当前小波系数对应的全局阈值λ_i：

其中median(ω_i)代表小波系数ω_i对应的中值，m₁,m₂代表图像的行和列，k代表系数；

(3)通过每一层小波系数ω_i和对应的全局阈值λ_i，求出每一层小波系数对应的估计小波系数；

其中，α＞1,r是调整因子，sgn(x)代表信号函数，当ω_i大于0时，信号函数值为1；当ω_i小于0，信号函数值为-1；

(4)使用小波重构算法得到去噪后的图像：针对每一层估计小波系数得到去噪后的图像，

其中N代表最大小波分解层数，M_N(x,y)代表经过小波分解算法后得到的低频分量，代表每一层估计小波系数。

进一步地，所述步骤(1)具体为：采用à trous小波和它的β-3 spline version对图像进行小波分解，包括如下步骤：

a.初始化i＝0，则第0层的图像为原图M₀；

b.变量i自加，图像M_i-1每行每列都与一个一维的核h进行卷积，卷积后图像表示为M_i,核h表示为矩阵并且在矩阵的元素之间插入(2^i-1-1)个零；

c.计算每一层小波系数：ω_i(k)＝M_i-1(k)-M_i(k)。

本发明的有益效果是：本发明具有如下优点：本发明提出的基于à trous小波阈值去噪算法，是对基于硬阈值的各向同性非抽取离散小波算法的改进，具体操作是使用特定的小波对高通量基因测序的图像进行小波分解，针对每一层小波系数，提出使用范数计算全局阈值。通过全局阈值和每一层小波系数构建估计小波系数表达式，最后使用小波重构算法得到去噪后的图像。本发明的方法满足Nyquist采样定理，具有移不变性；鲁棒性高，针对不同类型的噪声，算法的去噪效果都很显著；对比各种基于小波阈值去噪算法，图像得到的信噪比结果是最优的。

附图说明

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

图2是原始的高通量基因测序图像；

图3是添加标准差δ_n为20的高斯噪声的测序图；

图4是添加标准差δ_n为30的高斯噪声的测序图；

图5是添加泊松噪声后的测序图像，其中图5(a)是含噪声图像，图5(b)是图5(a)中的某亮斑灰度分布图，图5(c)是图5(a)中某亮斑去背景的灰度分布图；

图6是使用IUWT(hard thresholding)算法进行图像去噪的示意图，其中图6(a)是IUWT(hard thresholding)去噪图，图6(b)是图6(a)中的某亮斑灰度分布图，图6(c)是图6(a)中某亮斑去背景的灰度分布图；

图7是使用IUWT(soft thresholding)算法进行图像去噪的示意图，其中图7(a)是IUWT(soft thresholding)去噪图，图7(b)是图7(a)中的某亮斑灰度分布图，图7(c)是图7(a)中某亮斑去背景的灰度分布图；

图8是使用本发明的方法进行图像去噪的示意图，其中图8(a)是本发明的方法的去噪图，图8(b)是图8(a)中的某亮斑灰度分布图，图8(c)是图8(a)中某亮斑去背景的灰度分布图。

具体实施方式

下面结合附图说明及具体实施方式对本发明进一步说明。

本发明的主要步骤包括小波分解，阈值选取，估计去噪后的小波系数，小波重构等步骤，如附图1所示，本发明的针对高通量基因测序图像的小波去噪具体实现的过程如下：

(1)对测序图像使用特定的小波函数进行小波分解，得到每一层小波系数；

传统的小波基函数有许多种，比如离散小波变换，Mallet变换，和à trous小波变换，由于à trous小波满足Nyquist定理，具有移不变性，本发明优先选择à trous小波和它的β-3spline version对图像进行小波分解，其中可包括如下步骤：

1.1初始化i＝0，则第0层的图像为原图M₀；

1.2变量i自加，图像M_i-1每行每列都与一个一维的核h进行卷积，卷积后图像表示为M_i。核h表示为矩阵并且在矩阵的元素之间插入(2^i-1-1)个零；

1.3计算每一层小波系数：

ω_i(k)＝M_i-1(k)-M_i(k) (5)

(2)针对每一层小波系数，计算出当前小波系数对应的全局阈值；

通过步骤1得到每层小波系数ω_i，则计算得到该层对应的全局阈值：

式中median(ω_i)代表小波系数ω_i对应的中值，m₁,m₂代表图像的行和列，k代表系数。根据步骤(1)得到每一层小波系数，每一层小波系数的大小与真实图像的大小是相同的。由于使用范数计算出来的阈值会随噪声的影响而变化幅度较大，相比于范数，许多文献提出使用范数提出鲁棒性的特征提取算法。因此针对每一层小波系数，设计使用范数构造对应的阈值表达式。

(3)通过每一层小波系数和对应的全局阈值，求出每一层小波系数对应的估计小波系数；

式中α＞1,r是调整因子，sgn(x)代表信号函数，当ω_i大于0时，信号函数值为1；当ω_i小于0，信号函数值为-1。当小波系数ω_i的数值越大时，估计小波系数与真实小波系数之间数值越接近，说明当前小波系数代表真实信号时，估计小波系数会保留更多的真实信号的能量。反之，当小波系数ω_i的数值越小时，估计小波系数与真实小波系数之间差值越大，说明当前小波系数代表噪声信号时，估计小波系数会抑制噪声信号的能量。因此，本发明提出的估计小波系数表达式避免传统的小波软阈值去噪算法的缺点。

(4)使用小波重构算法得到去噪后的图像。

针对每一层估计小波系数得到去噪后的图像：

式中N代表最大小波分解层数，M_N(x,y)代表经过小波分解算法后得到的低频分量，代表每一层估计小波系数。

本发明将通过主观和客观两个方面验证算法的合理性和有效性。为了衡量各个算法去噪的效果，通过以下公式计算信噪比、均方误差作为评判标准：

上式中f(i,j)代表原图，f'(i,j)代表去噪后的图像。M,N代表图像的行高和列高。

首先在附图3、附图4、附图5(a)的测试图像中，使用本发明提供的方法、基于硬阈值的小波阈值去噪算法、基于软阈值小波阈值去噪算法三个算法分别进行去噪处理。各个算法计算得到的信噪比、均方误差结果分别存入表1、表2中。

表1不同算法在含噪声图像的SNR比较

表2不同算法在含噪声图像的MSE比较

表1和表2的实验结果表明，通过选取合适的参数，本发明的实验结果能够较好的提高图像的信噪比和降低图片的均方误差。在含有不同强度的噪声的图像中，本发明提出的模型计算得到的信噪比和均方误差的结果是最优的。本发明比IUWT(based on hardthresholding)算法得到的信噪比结果平均要高2dB左右，均方误差平均要小30左右；比IUWT(based on soft thresholding)算法得到的信噪比结果平均要高5dB左右，均方误差平均要小52左右。通过以上实验数据，说明本发明比传统的基于小波阈值去噪算法更优越，去噪性能更好。因此，本发明在降低图像高斯噪声，提高图像质量上具有合理性。

为了能从视觉上证明本发明的效果，通过对附图5(a)含有泊松噪声的图像进行实验，画出各种算法去噪后的图像。为了能够清楚的描绘去噪算法对细节特征的处理，对同一个区域的亮斑，画出不同算法去噪后的三维灰度曲面图像。由于使用基于à trous和它的β-3spline version小波对图像进行小波分解与重构处理，该小波是针对每一层小波系数进行两次平滑操作，将噪声图像的背景像素进行平滑。通过观察图6(b)、附图7(b)和附图8(b)发现，这三种算法都能够将含噪声图像中背景噪声像素统一，消除背景像素中泊松噪声，为以后碱基的识别提供了帮助。通过观察图6(c)、附图7(c)和附图8(c)，发现IUWT(softthresholding)的亮斑灰度值与原始的图像相比，重构后的灰度图与原图差别很大。而IUWT(hard thresholding)算法与本发明提出的模型重构的碱基图像的灰度值保留住大部分待测碱基的灰度信息，防止重构后的图像出现失真。并且本发明提出的模型使用范数计算得到阈值，而另外两种去噪算法使用范数计算得到阈值，对比各个算法降噪后的效果图(图6(c)、附图7(c)和附图8(c))说明本发明得到的去噪效果更具有鲁棒性和有效性，因为本发明重构的灰度图像能够较好的保留原始碱基图像中灰度信息，防止重构图像失真。

本发明的主要贡献在于：⑴使用à trous小波和它的β-3 spline version对图像进行小波分解与重构，算法具有移不变性；(2)相比于范数计算出来的阈值会随噪声的影响而变化幅度较大，使用范数计算每一层小波系数的全局阈值，提高算法的鲁棒性和去噪效果；(3)本发明花费的时间比基于局部阈值小波去噪算法要短，算法效率高；(4)克服传统的小波阈值去噪算法失真等缺点，本发明提出一种新的估计小波系数表达式,去噪后的图片效果更显著。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (1)

1.一种针对高通量基因测序图像的去噪方法，其包括以下步骤：

(1)对测序图像使用小波函数进行小波分解，得到每一层小波系数ω_i，具体为采用àtrous小波和它的β-3spline version对图像进行小波分解，包括如下步骤：

a.初始化i＝0，则第0层的图像为原图M₀；

b.变量i自加，图像M_i-1每行每列都与一个一维的核h进行卷积，卷积后图像表示为M_i，核h表示为矩阵并且在矩阵的元素之间插入(2^i-1-1)个零；

c.计算每一层小波系数：ω_i(k)＝M_i-1(k)-M_i(k)；

(2)针对每一层小波系数，计算出当前小波系数对应的全局阈值λ_i：

<mrow> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mo>=</mo> <mfrac> <mi>k</mi> <mrow> <msub> <mi>m</mi> <mn>1</mn> </msub> <mo>*</mo> <msub> <mi>m</mi> <mn>2</mn> </msub> </mrow> </mfrac> <mo>|</mo> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>m</mi> <mi>e</mi> <mi>d</mi> <mi>i</mi> <mi>a</mi> <mi>n</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msub> <mo>|</mo> <mn>1</mn> </msub> <mo>,</mo> </mrow>

其中median(ω_i)代表小波系数ω_i对应的中值，m₁,m₂代表图像的行和列，k代表系数；

(3)通过每一层小波系数ω_i和对应的全局阈值λ_i，求出每一层小波系数对应的估计小波系数；

<mrow> <msub> <mover> <mi>&amp;omega;</mi> <mo>~</mo> </mover> <mi>i</mi> </msub> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>0</mn> <mo>,</mo> <mrow> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>|</mo> </mrow> <mo>&amp;le;</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>sgn</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>-</mo> <mi>r</mi> <mo>*</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mo>*</mo> <msup> <mi>&amp;alpha;</mi> <mfrac> <mn>1</mn> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> <mo>|</mo> <msub> <mo>|</mo> <mn>2</mn> </msub> </mrow> </mfrac> </msup> <mo>)</mo> </mrow> <mo>,</mo> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>&gt;</mo> <msub> <mi>&amp;lambda;</mi> <mi>i</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow>

其中，α＞1,r是调整因子，sgn(x)代表信号函数，当ω_i大于0时，信号函数值为1；当ω_i小于0，信号函数值为-1；

(4)使用小波重构算法得到去噪后的图像：针对每一层估计小波系数得到去噪后的图像，

<mrow> <mi>M</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>M</mi> <mi>N</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>+</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msub> <mover> <mi>&amp;omega;</mi> <mo>~</mo> </mover> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>,</mo> </mrow>

其中N代表最大小波分解层数，M_N(x,y)代表经过小波分解算法后得到的低频分量，代表每一层估计小波系数。