一種基于偏微分方程的信號(hào)增強(qiáng)模型

李 俊 李遠(yuǎn)祿，2 蔣 民

（1.南京信息工程大學(xué)自動(dòng)化學(xué)院，南京，210044；2.江蘇省大氣環(huán)境與裝備技術(shù)協(xié)同創(chuàng)新中心，南京，210044）

引 言

峰檢測(cè)被應(yīng)用于很多領(lǐng)域，如心電圖、基于質(zhì)譜分析的癌癥診斷、水質(zhì)分析以及化學(xué)混合物鑒定等[1-2]。由于實(shí)驗(yàn)環(huán)境和儀器分辨水平等因素的影響，以及譜峰自身的特點(diǎn)，會(huì)存在部分弱峰，或者部分重疊峰，若直接對(duì)信號(hào)進(jìn)行峰檢測(cè)，容易造成峰的漏檢和錯(cuò)檢。因此，如何提高重疊峰的分辨力，或者對(duì)弱峰進(jìn)行增強(qiáng)是一個(gè)關(guān)鍵問(wèn)題[3]。

常見(jiàn)的信號(hào)增強(qiáng)方法包括Fourier去卷積方法、連續(xù)小波方法和導(dǎo)數(shù)譜方法。Kauppinen等人提出了Fourier去卷積法[4]；Du等人提出了連續(xù)小波的峰檢測(cè)方法[5]；導(dǎo)數(shù)譜方法是增強(qiáng)信號(hào)的常用方法[6]，目前被廣泛應(yīng)用于水體分析、乳腺腫瘤檢測(cè)和心電圖檢測(cè)等[7-9]。然而導(dǎo)數(shù)譜方法在增強(qiáng)譜峰的同時(shí)也會(huì)降低信號(hào)的信噪比，因此通常將導(dǎo)數(shù)譜方法和平滑方法相結(jié)合使用，如與小波方法結(jié)合用于化學(xué)信號(hào)檢測(cè)[10]，與Savitzky-Golay方法結(jié)合測(cè)定紅外光譜[11]，與高斯平滑方法結(jié)合用于蛋白質(zhì)譜峰檢測(cè)[12]等。

平滑方法較多，主要包括時(shí)域和頻域兩類方法。在時(shí)域中，最簡(jiǎn)單的平滑方法是滑動(dòng)均值濾波，即將相鄰的奇數(shù)個(gè)點(diǎn)求平均值代替原中心點(diǎn)[13]；Savitzky-Golay濾波目前被廣泛使用，該方法是一種廣義的滑動(dòng)均值濾波方法，它將一小組連續(xù)數(shù)據(jù)點(diǎn)做最小二乘擬合，并將多項(xiàng)式擬合曲線的中心點(diǎn)作為輸出[14],相比于滑動(dòng)均值濾波法，該方法具有更好的保峰效果；高斯濾波是對(duì)滑動(dòng)均值濾波平滑窗口的改進(jìn)方法，使用高斯函數(shù)作為平滑窗口[15]；類似地，還有使用kaiser窗作為平滑窗口的濾波方法[16]。在頻域中，需要先將信號(hào)轉(zhuǎn)換到頻域，再將高頻系數(shù)拋棄，然后重構(gòu)得到平滑信號(hào)。小波方法是目前最常用的頻域平滑方法，該方法可以較好地保護(hù)信號(hào)的特征，然而結(jié)果與所用小波基和尺度有關(guān)，不同尺度下所得結(jié)果可能相差很大[17-18]。

20世紀(jì)80年代，Witkin發(fā)現(xiàn)齊次線性擴(kuò)散方程的解等價(jià)于一定尺度下初始信號(hào)與高斯函數(shù)的卷積，因此偏微分方程在信號(hào)處理中越來(lái)越受到重視[19]。然而線性擴(kuò)散濾波在處理信號(hào)時(shí)沒(méi)有考慮信號(hào)的特征，在去除噪聲同時(shí)也可能會(huì)模糊信號(hào)的細(xì)節(jié)。針對(duì)這個(gè)問(wèn)題，Perona和Malik在1990年提出了經(jīng)典非線性擴(kuò)散方法，它根據(jù)原始信號(hào)的特征在信號(hào)的不同位置設(shè)置不同的擴(kuò)散強(qiáng)度，既有好的平滑效果，又能保護(hù)信號(hào)細(xì)節(jié)特征[20]。近年來(lái)，經(jīng)典非線性擴(kuò)散濾波被廣泛運(yùn)用到圖像處理中，在保護(hù)圖像紋理特征方面有很好的作用[21-23]。

鑒于非線性擴(kuò)散模型有保護(hù)信號(hào)細(xì)節(jié)特征的能力，本文將非線性擴(kuò)散與導(dǎo)數(shù)譜方法結(jié)合，得到一種新的信號(hào)增強(qiáng)模型。具體步驟是先對(duì)原信號(hào)應(yīng)用導(dǎo)數(shù)譜方法增強(qiáng)，再將增強(qiáng)后的信號(hào)作為非線性擴(kuò)散模型的初始信號(hào)，經(jīng)擴(kuò)散后得到平滑的增強(qiáng)信號(hào)。

1 基于非線性擴(kuò)散的導(dǎo)數(shù)譜增強(qiáng)模型

經(jīng)典非線性擴(kuò)散模型如下[20]

式中：f(x)為原始信號(hào)；g[u(x,t)]為擴(kuò)散函數(shù)。

導(dǎo)數(shù)譜增強(qiáng)模型如下[6]

式中：f為原始信號(hào)；F為增強(qiáng)后的信號(hào)；n為導(dǎo)數(shù)階次，通常取偶數(shù)；c為增強(qiáng)系數(shù)。

將經(jīng)典非線性擴(kuò)散模型與導(dǎo)數(shù)譜增強(qiáng)模型相結(jié)合，得到本文模型為

式中

目前，已有多種非線性擴(kuò)散方程的求解方法[21-22]，由于有限差分易于處理，并且實(shí)際數(shù)字信號(hào)已經(jīng)離散，因此本文采用有限差分格式求解。

有限差分格式主要有兩種方案，分別是顯式差分格式和隱式差分格式，由于隱式差分格式是無(wú)條件穩(wěn)定的，因此本文選取隱式差分格式對(duì)模型進(jìn)行離散。

隱式差分格式為

式中：i表示序號(hào)；k表示第k次迭代；τ為時(shí)間離散步長(zhǎng)；h為空間離散步長(zhǎng)。

將式（6）寫(xiě)成矩陣與向量的形式，有

式中

并且Bk-1是可逆的。

具體算法為：

（1）給定增強(qiáng)系數(shù)c，擴(kuò)散強(qiáng)度控制參數(shù)λ，迭代次數(shù)N，時(shí)間步長(zhǎng)τ，其中h可以取1；

（2）按式（2）得到原始信號(hào)的增強(qiáng)信號(hào)U0；

（3）按式（4）得到離散的gk-1i并計(jì)算βk-1i；

（4）構(gòu)造矩陣Bk-1；

（5）按式（7）迭代，即可得到平滑的增強(qiáng)信號(hào)。

2 模擬數(shù)據(jù)及結(jié)果討論

2.1 模擬數(shù)據(jù)

本文用洛倫茲峰來(lái)驗(yàn)證所提出的模型，它的產(chǎn)生方法如下

式中：n為峰的數(shù)目；Ai，μi和σi分別為第i個(gè)峰的峰高、峰位置和峰寬參數(shù)。

2.2 擴(kuò)散函數(shù)對(duì)擴(kuò)散過(guò)程的影響

非線性擴(kuò)散模型與線性擴(kuò)散模型的區(qū)別在于存在控制擴(kuò)散強(qiáng)度的擴(kuò)散函數(shù),如果擴(kuò)散函數(shù)為常數(shù)，則退化成線性擴(kuò)散模型，也就是常見(jiàn)的高斯平滑。它的缺點(diǎn)是在去除噪聲的同時(shí)會(huì)鈍化峰，尤其在信噪比較低時(shí)，容易形成重疊峰。非線性擴(kuò)散模型使用高斯函數(shù)作為擴(kuò)散函數(shù)，擴(kuò)散強(qiáng)度隨峰高的增加而減弱，從而到達(dá)保峰的目的。

在擴(kuò)散函數(shù)中，λ是控制擴(kuò)散強(qiáng)度的閾值。如果λ取值很大，擴(kuò)散函數(shù)的值趨近于1，這時(shí)相當(dāng)于高斯平滑。當(dāng)選擇合適的λ時(shí)，既可以保峰又可降噪。通常λ的值可根據(jù)經(jīng)驗(yàn)給定，當(dāng)然也可根據(jù)公式給定，若期待最高峰處的擴(kuò)散強(qiáng)度為一個(gè)小值ε，那么第k次迭代的λ可設(shè)為

當(dāng)噪聲水平較高時(shí)，需要在前幾次迭代中取較大的ε值，通常情況下取0.9以上；之后的迭代中則可取較小的ε值，從而較好地保護(hù)峰的形狀。

圖1是通過(guò)式(8)得到的由兩個(gè)不同高度的洛倫茲峰組成的模擬信號(hào)。對(duì)圖1中的噪聲信號(hào)做平滑，設(shè)較低峰處的擴(kuò)散強(qiáng)度為0.1，則較高峰處的擴(kuò)散強(qiáng)度ε為0，其擴(kuò)散函數(shù)如圖2(a)所示，平滑結(jié)果如圖2(b)所示，平滑后較高峰處仍存在較大噪聲；設(shè)較高峰處的擴(kuò)散強(qiáng)度ε為0.1，其擴(kuò)散函數(shù)如圖2(c)所示，平滑結(jié)果如圖2(d)所示，平滑后效果較好。由圖2可以看出，對(duì)最高峰取不同的擴(kuò)散強(qiáng)度ε，平滑效果是有差異的。

圖1 模擬信號(hào)Fig.1 Simulated signal

圖2 不同情況下擴(kuò)散函數(shù)對(duì)平滑結(jié)果的影響Fig.2 Influence of diffusion function on smoothing results under different conditions

2.3 信號(hào)產(chǎn)生與實(shí)驗(yàn)對(duì)比

為了對(duì)比增強(qiáng)與未增強(qiáng)的效果，本文通過(guò)式(8)的洛倫茲峰模型產(chǎn)生模擬信號(hào)，如圖3所示。圖3(a，b)分別是由8個(gè)和12個(gè)洛倫茲峰組成的信號(hào)，信號(hào)中存在弱峰和重疊峰，可以看出噪聲信號(hào)增強(qiáng)后噪聲水平也會(huì)被增強(qiáng)。

圖4是將圖3中2個(gè)信號(hào)經(jīng)過(guò)本文模型增強(qiáng)后得到的結(jié)果。為對(duì)比導(dǎo)數(shù)譜對(duì)峰的增強(qiáng)效果，需要為增強(qiáng)系數(shù)取一個(gè)合適值，本文對(duì)增強(qiáng)系數(shù)分別取0（未增強(qiáng)）和其他值時(shí)的增強(qiáng)效果進(jìn)行對(duì)比。實(shí)驗(yàn)中兩個(gè)信號(hào)取相同的平滑參數(shù)，時(shí)間步長(zhǎng)τ取1，迭代次數(shù)取20，前5次最高峰擴(kuò)散強(qiáng)度ε取0.99，之后取0.05。由圖4可以看出本文模型具有很好的峰增強(qiáng)效果，能增強(qiáng)弱峰的幅度和提高重疊峰的分離度。

為更好地表現(xiàn)增強(qiáng)的效果，本文將圖4的檢測(cè)結(jié)果用表1表示。通過(guò)圖4和表1的對(duì)比可以看出，只對(duì)噪聲信號(hào)做平滑，其中的弱峰會(huì)和噪聲容易一起被平滑掉，還容易導(dǎo)致重疊峰，從而影響檢測(cè)結(jié)果；使用本文模型處理后，信號(hào)中的弱峰幅度被增大，如信號(hào)1中x=150處，x=310處，信號(hào)2中x=270處，x=420處的峰，同時(shí)重疊峰的分辨會(huì)提高，如信號(hào)1中x=270處，信號(hào)2中x=160處，x=330處，x=380處的峰。

圖3 重疊峰真實(shí)信號(hào)、噪聲信號(hào)和增強(qiáng)后的噪聲信號(hào)（SNR=20 dB）Fig.3 Real signal,noisy signal and enhanced noisy signal of overlapping peaks(SNR=20 dB)

圖4 基于非線性擴(kuò)散的導(dǎo)數(shù)譜增強(qiáng)模型效果對(duì)比Fig.4 Comparison for the derivative spectrum enhancement models based on nonlinear diffusion

2.4 基于非線性擴(kuò)散的導(dǎo)數(shù)譜增強(qiáng)模型與其他模型的對(duì)比

導(dǎo)數(shù)譜方法常和小波方法（WMD）或Savitzky-Golay方法（SGD）結(jié)合應(yīng)用于信號(hào)增強(qiáng)，因此為進(jìn)一步驗(yàn)證基于非線性擴(kuò)散的導(dǎo)數(shù)譜增強(qiáng)模型（CDD）的效果，將其與前兩種模型作對(duì)比。在選取參數(shù)時(shí)盡量保證所有子峰都被檢測(cè)到，同時(shí)觀察信號(hào)中其他位置的平滑效果，兩個(gè)信號(hào)增強(qiáng)系數(shù)分別取100和150，WMD尺度取3，小波基選擇“sym12”，SGD框長(zhǎng)度取13，階次取2，迭代次數(shù)取10，3種模型對(duì)圖2兩個(gè)信號(hào)的增強(qiáng)結(jié)果如圖5所示。

通過(guò)兩組結(jié)果的對(duì)比，可以發(fā)現(xiàn)3種模型都具有較好的峰增強(qiáng)效果。然而和真實(shí)信號(hào)相比，WMD模型和SGD模型的增強(qiáng)結(jié)果中新增加了一些偽峰（與真實(shí)峰相比新產(chǎn)生的峰），由于這些偽峰的幅度較大，容易導(dǎo)致誤檢，如信號(hào)1中x=200處，信號(hào)2中x=200處和x=250處等。相比于其他兩種模型，本文模型在增強(qiáng)峰的同時(shí)，對(duì)峰以外的部分具有更好的去噪效果，有效降低了峰檢測(cè)中偽峰的數(shù)目。

2.5 基于非線性擴(kuò)散的導(dǎo)數(shù)譜增強(qiáng)模型的應(yīng)用

基于非線性擴(kuò)散的導(dǎo)數(shù)譜增強(qiáng)模型可以應(yīng)用于MALDI質(zhì)譜，MALDI質(zhì)譜的檢測(cè)對(duì)于發(fā)現(xiàn)疾病生物標(biāo)志物和研究藥物代謝等方面具有重要意義。本文從Bioinformatics Toolbox中選取一組MALDI模擬數(shù)據(jù)作為實(shí)驗(yàn)樣本(http：//bioinformatics.mdanderson.org/Supplements/Datasets/Simulations/index.html)，使用本文模型對(duì)其做處理。為方便對(duì)比，分別取增強(qiáng)系數(shù)=0（未增強(qiáng)）和增強(qiáng)系數(shù)=200，處理結(jié)果如圖6所示，其中豎線表示真實(shí)峰的位置。

從圖6可以看出本文模型對(duì)質(zhì)譜具有一定的增強(qiáng)效果，弱峰的高度被增加，同時(shí)從8 200左右的位置可以看出該模型能有效增加重疊峰的分離度。

表1 使用基于非線性擴(kuò)散的導(dǎo)數(shù)譜增強(qiáng)模型后峰檢測(cè)結(jié)果對(duì)比Tab.1 Comparison of peak detection performance using derivative spectrum enhancement model based on nonlinear diffusion

圖5 不同模型信號(hào)增強(qiáng)結(jié)果對(duì)比Fig.5 Comparison of enhancement results for different models

圖6 應(yīng)用增強(qiáng)與未增強(qiáng)質(zhì)譜的處理效果對(duì)比Fig.6 Comparison for the mass spectrometry with enhanced methods

3 結(jié)束語(yǔ)

本文提出了一種新的信號(hào)增強(qiáng)模型——基于非線性擴(kuò)散的導(dǎo)數(shù)譜增強(qiáng)模型，并通過(guò)模擬信號(hào)和MALDI質(zhì)譜驗(yàn)證所提模型的效果。結(jié)果表明所提模型能有效增強(qiáng)峰的特征，當(dāng)信號(hào)中存在噪聲，并有弱峰和重疊峰存在時(shí)，該模型能在去除噪聲的同時(shí)有效增強(qiáng)弱峰的幅度和提高重疊峰的分離度。