多峰曲線擬合

1、目錄1 背景介紹22 數(shù)據(jù)處理42.1 高斯擬合42.1.1 共軛梯度法52.1.2 0.618法62.1.3 計(jì)算效果72.2 分段多項(xiàng)式擬合82.2.1 多項(xiàng)式擬合的不足82.2.2 分段標(biāo)準(zhǔn)82.2.3 峰位、峰面積計(jì)算113 軟件介紹114 結(jié)語12多峰曲線擬合1 背景介紹工程上常采用電泳的方法來分離不同的DNA片段。帶電微粒在電場(chǎng)作用下，向著與其電性相反的電極移動(dòng)，稱為電泳(electrophoresis)。利用帶電粒子在電場(chǎng)中移動(dòng)速度不同而達(dá)到分離的技術(shù)稱為電泳技術(shù)。常用微流控芯片來做為DNA混合液的載體。微流控芯片技術(shù)是把生物、化學(xué)醫(yī)學(xué)分析過程的樣品制備、反應(yīng)、分離、檢測(cè)等基本操

2、作單元集成到一塊微米尺度的芯片上，自動(dòng)完成分析全過程。將DNA混合液放在特制的微流控芯片中，可以實(shí)現(xiàn)對(duì)溶液的精確控制，對(duì)DNA片段的檢測(cè)液變得方便簡(jiǎn)單。圖1為一片微流控芯片實(shí)物圖。實(shí)驗(yàn)室采用電泳技術(shù)來分離DNA片段。首先對(duì)溶液中的DNA片段進(jìn)行染色處理，然后將溶液注入微流控芯片。在微流控芯片的溝道兩端加載高壓進(jìn)行電泳，用激光照射經(jīng)染色處理的DNA片段使發(fā)出熒光，采用光電倍增管來檢測(cè)熒光的存在進(jìn)而檢測(cè)電泳過程,流程示意圖如圖2。用單片機(jī)采集光電倍增管的電流信號(hào)傳入計(jì)算機(jī)得到時(shí)間-電壓曲線（如圖3），根據(jù)曲線中峰的位置的個(gè)數(shù)來判斷溶液中DNA片段的種類。圖1 微流控芯片圖2 DNA檢測(cè)分析系統(tǒng)示意

3、圖圖3 電泳曲線在圖3中，每一個(gè)峰代表一個(gè)DNA片段，峰的位置和面積代表著DNA片段的種類信息。為了從曲線中獲得DNA信息，需要對(duì)曲線進(jìn)行相應(yīng)的處理。可以先對(duì)曲線作擬合處理，再由擬合得到的數(shù)學(xué)表達(dá)式來求得曲線中峰的位置和面積。本文討論兩種曲線擬合方式：1、高斯曲線擬合；2、多項(xiàng)式擬合。2 數(shù)據(jù)處理2.1 高斯擬合由于曲線中有多個(gè)峰，因此考慮用高斯函數(shù)和洛倫茲函數(shù)的線性組合來擬合曲線，預(yù)計(jì)會(huì)有好的結(jié)果1。記擬合函數(shù)為： (1)可以看到，由洛倫茲函數(shù)和高斯函數(shù)組合而成。其中各組合系數(shù)的意義如下：洛倫茲函數(shù)和高斯函數(shù)的組合比例，；：峰位置信息，；：峰的半高寬，；：峰強(qiáng)度信息，。要對(duì)電泳曲線進(jìn)行擬合，

4、就是要求參數(shù)，使函數(shù)能夠和曲線匹配的最好。本文采用下式來刻畫這種匹配效果： (2)其中為測(cè)得數(shù)據(jù)點(diǎn)。這樣，擬合就轉(zhuǎn)化為求的極小值問題： (3)這是一個(gè)求多元非線性函數(shù)的極值問題，目前沒有解析法直接求得，可以采用迭代方法逐步獲取該函數(shù)的極小值。本文采用共軛梯度法來求解。2.1.1 共軛梯度法共軛梯度法是最優(yōu)化中最常用的方法之一。2它具有算法簡(jiǎn)便、存儲(chǔ)需求小等優(yōu)點(diǎn)，十分適合于大規(guī)模優(yōu)化問題。共軛梯度法是介于最速下降法與牛頓法之間的一個(gè)方法，它僅需利用一階導(dǎo)數(shù)信息，但克服了最速下降法收斂慢的缺點(diǎn)，又避免了牛頓法需要存儲(chǔ)和計(jì)算Hesse矩陣并求逆的缺點(diǎn)。共軛梯度法不僅是解大型線性方程組最有用的方法之一

5、，也是解大型非線性最優(yōu)化問題最有效的算法之一。共軛梯度法最早是由Hestenes和Stiefel在1952年提出來的，用于解正定系數(shù)矩陣的線性方程組。在這個(gè)基礎(chǔ)上，F(xiàn)letcher和Reeves在1964年首先提出了解非線性最優(yōu)化問題的共軛梯度法。由于共軛梯度法不需要矩陣存儲(chǔ)，且有較快的收斂速度和二次終止性等優(yōu)點(diǎn)，現(xiàn)在共軛梯度法已經(jīng)廣泛地應(yīng)用于實(shí)際問題中。本文采用再開始共軛梯度法來作算法設(shè)計(jì)，確定共軛方向采用Fletcher-Reeves公式，一維搜索采用0.618法。以下簡(jiǎn)述算法設(shè)計(jì)步驟。記問題為，其中為n維向量，終止誤差為，為目標(biāo)向量。step1：初始步，給出初始點(diǎn)，終止誤差，。step2

6、：計(jì)算，如果，停止迭代，輸出；否則令。step3：線性搜索求，使得，并令，。step4：計(jì)算。若，令，轉(zhuǎn)step2；如果，停止迭代。step5：若，令，轉(zhuǎn)step2.step6：計(jì)算，。step7：若果，令，轉(zhuǎn)step2；否則轉(zhuǎn)step3。在第三步中，要求采用一維搜索獲取最佳步長(zhǎng)，本文采用0.618法來作精確線性搜索。2.1.2 0.618法0.618法是一種用于求單峰函數(shù)的極小值的分割法。其基本思想是通過取試探點(diǎn)和進(jìn)行函數(shù)值的比較，使包含極小點(diǎn)的搜索區(qū)間不斷縮短，當(dāng)區(qū)間長(zhǎng)度縮短到一定程度時(shí)，區(qū)間上各點(diǎn)的函數(shù)值均接近極小值，從而各點(diǎn)可以看作為極小點(diǎn)的近似。這種方法僅需計(jì)算函數(shù)值，不涉及導(dǎo)數(shù)，又

7、稱直接法。下面介紹0.618法的技術(shù)細(xì)節(jié)。設(shè)包含極小點(diǎn)的初始搜索區(qū)間為，設(shè)在上是凸函數(shù)。0.618法的基本思想是在搜索區(qū)間上選取兩個(gè)對(duì)稱點(diǎn)，其中是區(qū)間上的黃金分割點(diǎn)，是的對(duì)稱點(diǎn)，且，比較這兩點(diǎn)處的函數(shù)值和的大小，若，則說明極小值在區(qū)間內(nèi)，故取新區(qū)間為，若，則說明極小值在區(qū)間內(nèi)，故取新區(qū)間為。新區(qū)間的長(zhǎng)度為原區(qū)間長(zhǎng)度的0.618倍。新區(qū)間包含原區(qū)間中兩個(gè)對(duì)稱點(diǎn)中的一點(diǎn)，我們只要再選一個(gè)對(duì)稱點(diǎn)，并利用這兩個(gè)新對(duì)稱點(diǎn)處的函數(shù)值繼續(xù)比較，重復(fù)這個(gè)過程，最后即可確定極小點(diǎn)。下面給出具體算法步驟。step1：選取初始數(shù)據(jù)。確定初始搜索區(qū)間和精度要求。計(jì)算最初兩個(gè)試探點(diǎn)， 計(jì)算和，令。step2：比較目標(biāo)函

8、數(shù)值。若，轉(zhuǎn)step3，否則轉(zhuǎn)step4。step3：若，則停止計(jì)算，輸出；否則，令，計(jì)算，轉(zhuǎn)step2。step4：若，則停止計(jì)算，輸出，否則，令：計(jì)算，轉(zhuǎn)step2。2.1.3 計(jì)算效果本文采用C#作為編程語言來實(shí)現(xiàn)上述算法。當(dāng)進(jìn)行多次調(diào)試之后，發(fā)現(xiàn)共軛梯度法并不能解決問題。現(xiàn)總結(jié)原因如下。1、迭代初始值難以確定。共軛梯度法在開始之前要求設(shè)定初始迭代值。該方法對(duì)初始值的設(shè)定非常敏感，如果初始值設(shè)置的不夠好，這種方法的收斂將會(huì)變得非常慢，甚至不能收斂。由于每次進(jìn)行電泳的DNA片段有可能不同，峰值信息業(yè)會(huì)不斷變化，因此沒有一個(gè)好的方法來為算法提供初始值。2、在一維線性搜索方法中，首先要求給出搜

9、索區(qū)間。這個(gè)搜索區(qū)間也是由于問題的實(shí)際背景決定的。在本文中，不能確定峰值信息的取值區(qū)間，因此一維搜索只能采用純數(shù)學(xué)意義上的搜索區(qū)間，這樣搜索負(fù)擔(dān)將會(huì)變得非常沉重，迭代計(jì)算變得非常耗時(shí)。經(jīng)過再三權(quán)衡，筆者決定放棄這種方法?？赡芄P者對(duì)共軛梯度法的認(rèn)識(shí)還不夠，也可能這種方法不適合本問題的求解。因此筆者決定采用數(shù)值分析中的多項(xiàng)式擬合算法來獲取峰值信息。2.2 分段多項(xiàng)式擬合2.2.1 多項(xiàng)式擬合的不足多項(xiàng)式擬合算法簡(jiǎn)單，易于實(shí)現(xiàn)，是工程應(yīng)用中常用的擬合方法。多項(xiàng)式擬合對(duì)于單峰、變化平緩的曲線擬合效果較好，但對(duì)于多峰曲線擬合來說，則差強(qiáng)人意。圖4是對(duì)電泳數(shù)據(jù)作40階多項(xiàng)式擬合的結(jié)果。從圖4可以看到，多項(xiàng)

10、式函數(shù)不能表現(xiàn)曲線的突變趨勢(shì)，因而便不能直接用來擬合多峰曲線。圖4 多項(xiàng)式擬合電泳曲線為了改變這種現(xiàn)狀，本文采用分段多項(xiàng)式擬合技術(shù)，實(shí)踐證明，將數(shù)據(jù)人為分段之后，對(duì)每一段單獨(dú)進(jìn)行擬合，效果很好3。分段的目的就是要把曲線平坦的部分和有峰的部分區(qū)分開，這樣才會(huì)得到好的擬合效果。下面介紹針對(duì)于本問題的一種分段方式。2.2.2 分段標(biāo)準(zhǔn)對(duì)于整段曲線來說，有峰的地方所占的比例不大，因此整段曲線的平均值將接近于背景暗電流值。在經(jīng)過試驗(yàn)以后，發(fā)現(xiàn)曲線的平均值如圖5所示。在圖5中，紅線表示了整段曲線的平均值。可以看到，這條紅線并不能“經(jīng)過”所有的峰，這樣就不能夠獲取分組依據(jù)。圖5 曲線平均值示意圖考慮到去掉一

11、些暗電流值，再將余下的數(shù)據(jù)作平均，則可能使得新的平均值“經(jīng)過”所有的峰。本文采用這樣一種策略：以第一次求得整條曲線的平均值作為標(biāo)準(zhǔn)，將每個(gè)數(shù)據(jù)與之比較，去掉所有小于平均值的數(shù)據(jù)點(diǎn)，對(duì)剩下的數(shù)據(jù)點(diǎn)再做一次平均，將新獲得的平均值作為分組依據(jù)。圖6顯示了新平均值的位置。可以看到，新的平均值經(jīng)過所有峰，這樣就可以根據(jù)新的平均值與原曲線的“交點(diǎn)”來劃分曲線了。圖6 新平均值示意圖對(duì)每一段曲線作40階多項(xiàng)式插值，最后得到的插值效果如圖7所示。由于多項(xiàng)式插值的算法已經(jīng)非常成熟，在此就不作詳細(xì)介紹4。需要指明的是，本文在實(shí)現(xiàn)多項(xiàng)式插值算法的時(shí)候，求解法方程組所采用的算法是Gauss按比例列主元消去法。圖8、圖

12、9為圖7采取局部放大之后的圖像，來顯示多項(xiàng)式插值的細(xì)節(jié)。圖7 分段多項(xiàng)式擬合效果圖8 擬合細(xì)節(jié)圖9 擬合細(xì)節(jié)2.2.3 峰位、峰面積計(jì)算在分段多項(xiàng)式插值完成后，每一個(gè)峰都可以用一個(gè)多項(xiàng)式函數(shù)來表達(dá)，那么對(duì)峰位的求解就轉(zhuǎn)化為求這個(gè)多項(xiàng)式函數(shù)在閉區(qū)間內(nèi)的最小值，對(duì)峰面積的求解就轉(zhuǎn)化為對(duì)此多項(xiàng)式函數(shù)在閉區(qū)間內(nèi)求定積分，這在算法的實(shí)現(xiàn)上都非常簡(jiǎn)單，使得程序設(shè)計(jì)變得非常輕松。圖10下方文本框中顯示了峰位、峰面積計(jì)算結(jié)果。圖10 峰位、峰面積計(jì)算結(jié)果3 軟件介紹本軟件采用微軟的WPF程序設(shè)計(jì)框架設(shè)計(jì)而成。這使得界面設(shè)計(jì)變得靈活簡(jiǎn)單，界面元素也變得更加漂亮。圖像顯示控件采用一款開源圖表控件ZedGraph，

13、這款控件編程簡(jiǎn)單，性能優(yōu)良。整個(gè)軟件界面分成三個(gè)部分。最上面是圖像顯示部分，用來顯示曲線。左下角為操作部分，有兩個(gè)按鈕，第一個(gè)用來載入數(shù)據(jù)，數(shù)據(jù)文件格式為.txt格式，第二個(gè)按鈕用來執(zhí)行擬合操作。右下角有一個(gè)文本框，用來顯示計(jì)算結(jié)果。具體構(gòu)造可以參看圖10。4 結(jié)語本文為了獲取曲線的峰值信息，嘗試了多種不同的方式進(jìn)行求解。實(shí)踐證明，復(fù)雜的算法并不一定能夠產(chǎn)生好的效果。就像本文之前介紹的共軛梯度法，用這個(gè)方法可以解決大型非線性最優(yōu)化問題，當(dāng)然也應(yīng)該能夠解決高斯擬合參數(shù)求解問題，但是在實(shí)際操作中由于此種方法對(duì)于初始迭代值的要求過高，沒有能夠成功實(shí)現(xiàn)，并且該算法計(jì)算復(fù)雜，耗時(shí)很長(zhǎng)，對(duì)于一個(gè)擬合問題來說，的確有大材小用之嫌。所以，本文隨后采用了熟悉的方法：多項(xiàng)式擬合。只是在擬合之前對(duì)數(shù)據(jù)作了分段處理。這樣一來，可以很完美的解決問題，而算法設(shè)計(jì)又非常簡(jiǎn)單，且計(jì)算非常輕松。所以，在考慮問題解決方式的時(shí)候，對(duì)于同樣能解決問題的方法，應(yīng)當(dāng)選擇那