局部脑血流的测定摘要无论对病人还是对医院而言，脑血流的测定是一种非常重要的技术手段，测定脑血流量可为研究人脑在不同的病理和生理条件下的功能提供客观指标，它对研究脑循环药物的药理作用也很有帮助，因此对于脑血流的研究具有重要的意义。

问题1，首先根据题设可知：一方面，由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比，并且其比例系数反映了该处的脑血流量；另一方面，动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比，由上述两方面可得到头部记数率关于脑部血流系数和呼出气记数率的表达式；再分析试验得到的呼出气记数率的数据，用Matlab 的拟合工具箱拟合得到关于呼出气记数率的函数；最后代入上述表达式求解微分方程即可得到关于脑部血流系数的表达式，表达式为)(5.11000)(5.1Kt te e K k t N ----=。

问题2，先用Matlab 拟合工具箱求出头部记数率的函数，用对比系数法可得到K 和k 的近似值，即3977.0,5015.0≈≈k K ，但这种方法是不精确的，只是用于后面方法得到参数的验证；将模型视为含两个参数的一次函数，对X ，Y 对数值进行拟合得参数K =0.5006，k =0.3937，对得到的数据进行误差分析可知脑血流的预测值和实际值波动比较大，因此这种算法也不太精确；再利用最小二乘法求解拟合后曲线参数的函数，由问题1得到的关于脑部血流系数的表达式即可得到脑血流系数：0.5000=K ，4001.0=k ，最后对得到的值进行误差分析，可知脑血流的预测值和实际值很吻合，比较符合题意。

关键词：脑血流系数 最小二乘法 曲线拟合 Matlab一问题重述用放射性同位素测定大脑局部脑血流量的方法如下：由受试者吸入含有某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处的放射性记数率(简称记数率)，同时测量他呼出气的记数率。

由于动脉血将肺部的放射性同位素传送至大脑，使脑部同位素增加，而脑血流又将同位素带离，使同位素减少。

实验证明由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比，其比例系数反映了该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。

动脉血从肺输送同位素至大脑引起脑部记数率上升的速度与当时呼出气的记数率成正比。

若某受试者的测试数据如附表1所示：根据以上题目所给的条件及数据，回答以下问题：1、建立确定脑部血流系数的数学模型；2、计算上述受试者的脑血流系数。

二模型假设1.脑部记数率的上升只与从肺部输送的放射性同位素有关，并且上升的速率与当时呼出气的记数率成正比；2.脑部记数率的下降只与当时该处的脑血流量有关，并且下降的速率与当时该处的记数率成正比，这里忽略了放射性同位素自身的衰变和其它因素；3.脑血流量在测定期间恒定，心脏博动、被测试者大脑活动、情感波动等带来的变化忽略不予考虑；4.每次仪器测量为相互独立事件，各测量值无记忆相关；5.放射性同位素在人体内传递是从吸入气体(含有放射物)开始的，并假定一次吸入，因此认为同位素在肺中瞬时达到最大浓度；6.在吸入气体瞬时，脑中放射物记数率为零；7.脑血流量与脑血流量系数成单值函数关系，求得后者即可确定前者。

三符号说明四 问题分析问题1，首先根据题设可知：一方面，由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比, 并且其比例系数反映了该处的脑血流量；另一方面，动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比，由上述两方面可得到头部记数率关于脑部血流系数和呼出气记数率的表达式；再分析试验得到的呼出气记数率的数据，用MATLAB 的拟合工具箱拟合得到关于呼出气记数率的函数；最后代入上述表达式求解微分方程即可得到关于脑部血流系数的表达式。

问题2，利用Matlab 里专门求解拟合后曲线参数的函数，再由问题1得到的关于脑部血流系数的表达式即可得到脑血流系数。

最后再对得到的值进行误差分析。

五 模型的建立与求解5.1求解脑血流系数的数学模型设某时刻0≥t 时，脑部记数率为)(t N ，在t ∆时刻后记数率为)(t t N ∆+，由题设及基本假设1和2可知，脑部记数率的增量)()(t N t t N N -∆+=∆只与下面两个因素有关：⑴动脉血从肺部输送放射性同位素至大脑引起脑部记数率的增量为1N ∆； ⑵脑血流将放射性同位素带离使得脑部记数率的减量为2N ∆。

综上，由医学实验可得：)()(21t KN dtdN t kP dtdN ==， 又)()()(21t N t N t N ∆-∆=∆所以有：dtdN dt dN dt dN 21-= 即)()(t KN t kP dtdN-= (1) 分析式(1)，要确定脑血流系数的模型，必须分析)(t P 和)(t N 的实验数据，观察其变化趋势。

首先用Matlab 绘出)(t P 和)(t N 的散点图并观察其变化趋势，)(t P 和)(t N 的散点图如下：图5-1 呼出气记数率的散点图图5-2 脑部记数率的散点图由上图可知，t 和)(t P 有近似于btae 的关系，而t 和)(t N 的关系暂时不能直接观察出，设btae t P =)(，用MATLAB 的拟合工具箱可得到参数a =1000，=b -1.5，而其相关系数为1，说明拟合的非常精确。

拟合后的图像如下图所示：图5-2 拟合函数)(t P 的图像由基本假设6，即0)0(=N 和式(1)联立可得一带有初值的微分方程：⎪⎩⎪⎨⎧=-=-0)0()(10005.1N t KN ke dtdN t (2) 解此微分方程得到)(5.11000)(5.1Kt te e K k t N ----= (3)5.2 求解脑血流系数的算法模型以下为求解)(t N 的算法模型：算法模型I ：再分析)(t N 的试验数据，同样用Matlab 的拟合工具箱里现有的指数函数可精确拟合出)(t N 的函数，表达式为：)(t N =t t e e 479.15015.039834032--⨯-⨯，可近似认为5.1479.1,39834032≈≈，则对比系数可得到3977.0,5015.0≈≈k K ，但这种方法是不精确的，只是用于后面方法得到参数的验证。

对于参数K 和k 的求解，可以利用Matlab 里专门求解拟合后曲线参数的函数lsqcurvefit 求得，0.4001k 0.5000, ==K 。

最后对得到的数据进行误差分析，把求得的K 和k 的值代入（3）式即可得到)(t N 的函数，求出预测值，再和实际值进行作差，用Matlab 绘出作差后的值和时间t 的关系图，如下所示：图5-3 误差分析比较图由图可知实际值与预测值的差值很小，在0上下波动，比较符合题意。

算法模型II ：由（1）变形可得：k t P t N t P t +⎪⎪⎭⎫ ⎝⎛-⨯=)()(K )()(N' （4） 即：k +⨯=X K Y （5）由（5）可以看出，因为K,k 为常数，所以Y 与X 在数值上应满足线性关系， 利用Matlab 拟合工具箱得到N(t)随时间变化关系为：t t e e t 479.15015.039834032(N ---=）（6） -1.479t-0.5015t5890.857e-2022.048e(t)N' += （7）由Matlab 编程求解，得到X ，Y 序列值，作图如下：图5-4 X-Y 散点图用Matlab拟合工具按（5）线性拟合如下图：图5-5 X-Y线性拟合并得到（5）式中参数：K=5006.0k=.03937对得到的数据进行误差分析，把求得的K和k的值代入（3）式即可得到)(tN的函数，求出预测值，再和实际值进行作差，用Matlab绘出作差后的值和时间t 的关系图，如下所示：图5-6 误差分析比较图由图可知实际值与预测值的差值比较大，因此此模型不太合理。

六模型的评价与推广1在建模时忽略了同位数的衰变以及动脉血从肺部到脑部所需要的时间，如在模型考虑这些因素后，只须在测试中测得这些因素的数值，用上述方法仍是容易实现的。

2本题模型比较简单，利用解微分方程和最小二乘法求得的脑部血流系数比较精确，拟合出的)(tN和)P的函数精度很高，第I种算法模型比较稳定。

(t七参考文献[1]曹卫华,郭止.最优化技术方法及MATLAB的实现[M],北京：化学工业出版社，2005.1[2]王家文，王皓，刘海.MATLAB7.0编程基础[M],北京：机械工业出版社，2005.7[3]刘志平，石林英.最小二乘法原理及其MATLAB实现[J],中国科技西部，2008，17(7):33-34附录附录1附录2x=[1.00:0.25:10.0];y=[1534,1528,1468,1378,1272,1162,1052,947,848,757,674,599,531,471,417,369,326,2 88,255,225,199,175,155,137,121,107,94,83,73,65,57,50,44,39,35,31,27];a0=[0 0];a=lsqcurvefit(@(a,x)(10000*a(1)*(exp(-1.5*x)-exp(-a(2)*x)))/(a(2)-1.5),a0,x,y)附录3x=[1.00:0.25:10.0];for i=1:37y=10000*0.4001*(exp(-1.5*x(i))-exp(-a(2)*x(i)))/(0.5000-1.5)end附录4t=1:0.25:10;N=[1534 1528 1468 1378 1272 1162 1052 947 848 757 674 599 531 471 417 369 326 288 255 225 199 175 155 137 121 107 94 83 73 65 57 50 44 39 35 31 27];P=[2231 1534 1054 724 498 342 235 162 111 76 52 36 25 17 12 8 6 4 3 2 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0];for i=1:37Nff(i) =-2022.048*exp(-0.5015*t(i))+5890.857*exp(-1.479*t(i));endfor i=1:37y(i)=Nff(i)/P(i);x(i)=-N(i)/P(i);end。