多峰曲线拟合

目录1 背景介绍 22 数据处理 42.1 高斯拟合 42.1.1共轭梯度法 52.1.20.618 法 62.1.3计算效果 72.2 分段多项式拟合 82.2.1多项式拟合的不足 82.2.2分段标准 82.2.3峰位、峰面积计算 113 软件介绍 114 结语 12多峰曲线拟合1背景介绍工程上常采用电泳的方法来分离不同的DNA片段。带电微粒在电场作用下， 向着与其电性相反的电极移动，称为电泳(electrophoresis)。利用带电粒子在电场 中移动速度不同而达到分离的技术称为电泳技术。常用微流控芯片来做为 DNA 混合液的载体。微流控芯片技术是把生物、化 学医学分析过程的样品制备、反应、分离、检测等基本操作单元集成到一块微米 尺度的芯片上，自动完成分析全过程。将DNA混合液放在特制的微流控芯片中， 可以实现对溶液的精确控制，对DNA片段的检测液变得方便简单。图1为一片 微流控芯片实物图。实验室采用电泳技术来分离DNA片段。首先对溶液中的DNA片段进行染色 处理，然后将溶液注入微流控芯片。在微流控芯片的沟道两端加载高压进行电泳， 用激光照射经染色处理的 DNA 片段使发出荧光，采用光电倍增管来检测荧光的 存在进而检测电泳过程,流程示意图如图2。用单片机采集光电倍增管的电流信号 传入计算机得到时间-电压曲线(如图 3)，根据曲线中峰的位置的个数来判断溶 液中DNA片段的种类。图1 微流控芯片激光机 片 自-21.81.61.41.210.80.6x 104图 2 DNA 检测分析系统示意图0 500 1000 1500 2000 2500 3000time/ms在图3中，每一个峰代表一个DNA片段，峰的位置和面积代表着DNA片段的种类信息。为了从曲线中获得DNA信息，需要对曲线进行相应的处理。可以先对曲线作拟合处理，再由拟合得到的数学表达式来求得曲线中峰的位置和面积本文讨论两种曲线拟合方式：1、高斯曲线拟合；2、多项式拟合。2数据处理2.1 高斯拟合由于曲线中有多个峰，因此考虑用高斯函数和洛伦兹函数的线性组合来拟合 曲线，预计会有好的结果1。记拟合函数为：vj 、 x 一 v j wj丿1可以看到，F(x)由洛伦兹函数1 + x 2F (x) = j=1(1-c )j1 +x-v._(j )2+ c I e wjjj和高斯函数e - x2组合而成。其中各组合系数的意义如下：C :洛伦兹函数和高斯函数的组合比例，C 0,1； jjV ：峰位置信息，V (0, +8； jjw ：峰的半高宽， w (0,+8； jjI ：峰强度信息， I (0,+8。jj要对电泳曲线进行拟合，就是要求参数c，v，w，I使函数F(x)能够和 jjjjQ=mi=1 j =1曲线匹配的最好。本文采用下式来刻画这种匹配效果：y - F(x )2 i = 1,2,m j = 1,2,n icj ,vj,wj ,Ij i其中(x , y )为测得数据点。这样，拟合就转化为求Q的极小值问题:iimin Q (c , v , w , I ) j = 1,2,nj j j j这是一个求多元非线性函数的极值问题，目前没有解析法直接求得，可以采用迭代方法逐步获取该函数的极小值。本文采用共轭梯度法来求解。2.1.1 共轭梯度法共轭梯度法是最优化中最常用的方法之一。2它具有算法简便、存储需求小 等优点，十分适合于大规模优化问题。共轭梯度法是介于最速下降法与牛顿法之 间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点， 又避免了牛顿法需要存储和计算 Hesse 矩阵并求逆的缺点。共轭梯度法不仅是解 大型线性方程组最有用的方法之一，也是解大型非线性最优化问题最有效的算法 之一。共轭梯度法最早是由Hestenes和Stiefel在1952年提出来的，用于解正定 系数矩阵的线性方程组。在这个基础上，Fletcher和Reeves在1964年首先提出 了解非线性最优化问题的共轭梯度法。由于共轭梯度法不需要矩阵存储，且有较 快的收敛速度和二次终止性等优点，现在共轭梯度法已经广泛地应用于实际问题 中。本文采用再开始共轭梯度法来作算法设计，确定共轭方向采用 Fletcher-Reev es 公式，一维搜索采用 0.618 法。以下简述算法设计步骤。/ 、 *记问题为min f x，其中?为n维向量，终止误差为 0， g二f， ?为 k丿目标向量。step1:初始步，给出初始点X。，终止误差0，k = 0。step2 :计算g = g (x )，如果g II ，停止迭代，输出g* =g0 ；否则令000 II0d 0 = -gU 0step3 :线性搜索求九，k使得f (x +九d ) =mink k kkkx = x + X d ， k = k +1。 k +1k k kstep4：计算gk令X = X0k，转 step2；如果 |g | 0，令 x = x，转 step2；否则转 step3。k k0 k在第三步中，要求米用一维搜索获取最佳步长九，本文米用0.618法来作精 k确线性搜索。2.1.20.618法0.618 法是一种用于求单峰函数的极小值的分割法。其基本思想是通过取试 探点和进行函数值的比较，使包含极小点的搜索区间不断缩短，当区间长度缩短 到一定程度时，区间上各点的函数值均接近极小值，从而各点可以看作为极小点 的近似。这种方法仅需计算函数值，不涉及导数，又称直接法。下面介绍 0.618 法的技术细节。设包含极小点的初始搜索区间为a,b，设(a)= f (x +ad )在a,b上 kk是凸函数。0.618法的基本思想是在搜索区间a,b上选取两个对称点a，a，12其中a是区间a,b上的黄金分割点，a是a的对称点，且a 申(a)，则说明极小值在区间1 2 1 2a ,b内，故取新区间为a ,b，若申(a ) 0。计算最11初两个试探点九，卩，11九=a + 0.328(b a )1 1 1 1卩二 a + 0.618（b a ）1 1 1 1计算甲（九）和申（卩），令k = 1。11step2:比较目标函数值。若申（九） 申（卩），转step3，否则转step4。kkstep3：若b -X 5，则停止计算，输出卩；否则，令k kka = X ,b = b ,X =卩，k+1k k +1 k k +1k申（X ）=申（卩）,卩 =a + 0.618（b - a ）k +1k k +1k +1k +1k +1计算甲（卩），转step2。k+1step4：若卩-a 5，则停止计算，输出X，否则，令：k kka = a ,b =卩,卩=Xk +1k k +1k k +1k申（卩）=申（X ）, X = a + 0.382（b a ）k +1k k +1k +1k +1 k +1计算申（X ），转step2。k+12.1.3 计算效果本文采用C#作为编程语言来实现上述算法。当进行多次调试之后，发现共轭 梯度法并不能解决问题。现总结原因如下。1、迭代初始值难以确定。共轭梯度法在开始之前要求设定初始迭代值。该 方法对初始值的设定非常敏感，如果初始值设置的不够好，这种方法的收敛将会 变得非常慢，甚至不能收敛。由于每次进行电泳的 DNA 片段有可能不同，峰值 信息业会不断变化，因此没有一个好的方法来为算法提供初始值。2、在一维线性搜索方法中，首先要求给出搜索区间。这个搜索区间也是由 于问题的实际背景决定的。在本文中，不能确定峰值信息的取值区间，因此一维 搜索只能采用纯数学意义上的搜索区间，这样搜索负担将会变得非常沉重，迭代 计算变得非常耗时。经过再三权衡，笔者决定放弃这种方法。可能笔者对共轭梯度法的认识还不 够，也可能这种方法不适合本问题的求解。因此笔者决定采用数值分析中的 多项式拟合算法来获取峰值信息。2.2 分段多项式拟合2.2.1 多项式拟合的不足多项式拟合算法简单，易于实现，是工程应用中常用的拟合方法。多项式拟 合对于单峰、变化平缓的曲线拟合效果较好，但对于多峰曲线拟合来说，则差强 人意。图 4 是对电泳数据作 40 阶多项式拟合的结果。从图 4 可以看到，多项式 函数不能表现曲线的突变趋势，因而便不能直接用来拟合多峰曲线。标題I 一 廿生茁 1050010001500200025003M0timifi fEra图4多项式拟合电泳曲线为了改变这种现状，本文采用分段多项式拟合技术，实践证明，将数据人为 分段之后，对每一段单独进行拟合，效果很好3。分段的目的就是要把曲线平坦 的部分和有峰的部分区分开，这样才会得到好的拟合效果。下面介绍针对于本问 题的一种分段方式。2.2.2 分段标准对于整段曲线来说，有峰的地方所占的比例不大，因此整段曲线的平均值将 接近于背景暗电流值。在经过试验以后，发现曲线的平均值如图 5 所示。在图 5中，红线表示了整段曲线的平均值。可以看到，这条红线并不能“经过”所有的 峰，这样就不能够获取分组依据。datal标题4 4I111111111111111111111111111105001000150020002M03D0time图5曲线平均值示意图考虑到去掉一些暗电流值，再将余下的数据作平均，则可能使得新的平均值 “经过”所有的峰。本文采用这样一种策略：以第一次求得整条曲线的平均值作 为标准，将每个数据与之比较，去掉所有小于平均值的数据点，对剩下的数据点 再做一次平均，将新获得的平均值作为分组依据。图 6 显示了新平均值的位置。 可以看到，新的平均值经过所有峰，这样就可以根据新的平均值与原曲线的“交 点”来划分曲线了。标题22 -J111J11J1111!111!20 -图6 新平均值示意图对每一段曲线作40阶多项式插值，最后得到的插值效果如图7所示。由于多项式插值的算法已经非常成熟，在此就不作详细介绍4。需要指明的是，本文在实现多项式插值算法的时候，求解法方程组所采用的算法是Gauss按比例列主元消去法。图8、图9为图 7 采取局部放大之后的图像，来显示多项式插值的细 节。标题22 H111!11!1n111!1111AHI - n 1图 7 分段多项式拟合效果曷丁由I1-4001500图 8 拟合细节dotal10.610.4145015001550time /rns图 9 拟合细节2.2.3峰位、峰面积计算在分段多项式插值完成后，每一个峰都可以用一个多项式函数来表达，那么 对峰位的求解就转化为求这个多项式函数在闭区间内的最小值，对峰面积的求解 就转化为对此多项式函数在闭区间内求定积分，这在算法的实现上都非常简单， 使得程序设计变得非常轻松。图10下方文本框中显示了峰位、峰面积计算结果。01000150020002500标盏I 臼创富 伽I3000Eime /ms图10 峰位、峰面积计算结果3 软件介绍本软件采用微软的 WPF 程序设计框架设计而成。这使得界面设计变得灵活 简单，界面元素也变得更加漂亮。图像显示控件采用一款开源图表控件ZedGraph, 这款控件编程简单，性能优良。整个软件界面分成三个部分。最上面是图像显示 部分，用来显示曲线。左下角为操作部分，有两个按钮，第一个用来载入数据， 数据文件格式为.tx格式，第二个按钮用来执行拟合操作。右下角有一个文本框， 用来显示计算结果。具体构造可以参看图10。4 结语本文为了获取曲线的峰值信息，尝试了多种不同的方式进行求解。实践证明， 复杂的算法并不一定能够产生好的效果。就像本文之前介绍的共轭梯度法，用这 个方法可以解决大型非线性最优化问题，当然也应该能够解决高斯拟合参数求解 问题，但是在实际操作中由于此种方法对于初始迭代值的要求过高，没有能够成 功实现，并且该算法计算复杂，耗时很长，对于一个拟合问题来说，的确有大材 小用之嫌。所以，本文随后采用了熟悉的方法：多项式拟合。只是在拟合之前对 数据作了分段处理。这样一来，可以很完美的解决问题，而算法设计又非常简单， 且计算非常轻松。所以，在考虑问题解决方式的时候，对于同样能解决问题的方 法，应当选择那些简单而易于实现的方法，这样才能得到最快的完成效率。参考文献：1 徐怡庄. 红外光谱曲线分峰拟合程序的编写. 光谱学与光谱分析， 1997,2.2 戴彧虹，袁亚湘 . 非线性共轭梯度法. 上海科学技术出版社. 2000,10.3 姜云，康智慧，王丛石，高锦岳. 用最小二乘法对多峰值数据进行曲线拟合及 曲线峰中心位置的精确求解. 物理实验. 2003.23.4 林成森 . 数值分析 . 科学出版社 . 2007.1.