局部脑血流的测定一. 问题简介脑血流量是诊断和治疗脑梗塞,脑出血,动脉瘤和先天性动,静脉血管畸形等脑血管疾病的主要依据。

测定脑血流量可为研究人脑在不同的病理和生理条件下的功能提供客观指标，它对研究脑循环药物的药理作用也很有帮助。

所以人们长期致力于寻找有效地测定脑血流量的方法。

近年来出现了以放射性同位素作示踪剂测定人脑局部血流量的方法。

这种方法大致可描述如下:由受试者吸入某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处放射性同位素的计数率(简称计数率)，同时测量他呼出气的计数率。

由于动脉血将肺部的放射性同位素输送至大脑，使脑部同位素增加，而脑血流又将同位素带离，使同位素减少，实验证明由脑血流引起局部地区计数率下降的速率与当时该处的记数率成正比，其比例系数反映了该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。

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

若某受试者的测试数据如下:时间(分) 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75头部记数率1534 1528 1468 1378 1272 1162 1052 947 848 757 674 599呼出气记数率2231 1534 1054 724 498 342 235 162 111 76 52 36时间(分) 4.00 4.25 4.50 4.75 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75头部记数率531 471 417 369 326 288 255 225 199 175 155 137呼出气记数率25 17 12 8 6 4 3 2 1 1 1 1时间(分) 7.00 7.25 7.50 7.75 8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.0头部记数率121 107 94 83 73 65 57 50 44 39 35 31 27呼出气记数率0 0 0 0 0 0 0 0 0 0 0 0 0试建立确定脑血流系数的数学模型并计算上述受试者的脑血流系数。

备注:该题目是上海市(1990 年)大学生数学建模竞赛A题。

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

三. 模型的建立与分析由于已知脑局部同位素的增减与已测定的头部记数率ht()和呼出气记数率p()t 成正比关系,于是很自然地确定以脑部同位素量,即脑部记数率作为讨论对象.。

1.原始模型的建立设某时刻t ≥0时,头部记数率为ht(),在Δt 时段后记数率ht()+Δt ，由假定可知, 头部记数率的增量Δhht= ()(+Δt -ht)仅与三个因素有关：(i) 肺动脉血将肺部的放射性同位素送至大脑,使脑部记数率增量为Δh ；(ii) 脑血流将同位素带离,脑记数率下降为Δh ；(iii)放射性同位素自身有衰减引起记数率下降量为Δh ,设其半衰期为τ.因此,由医学试验及假定有dh dh dh ln21 2 3= kp(),t = Kh(),t =- ht(),dt dt dt τ而Δht() = Δh ()t -Δh ()t +Δh ()t ，123于是dh dh dh dh123=-+ ,dt dt dt dt即dh ln2=-Kht()+kpt()- ht(). (dt τ133由于在测试时放射性同位素(如Xe)的半衰期τ一般很大,而测试时间又很短(大约十几分钟左右),由此假定τ→+∞,于是(1)式变为:dh=-Kh()t +kp()t . (dt2. 算法模型的建立与改进在建立算法模型之前,首先必须对p()t 进行预测。

作p()tt~ 和-λt ln pt( ) ~ t 的离散图(图1和图2),由此发现p()t 与t 有近似于Ae 的函数关系。

通过对ln pt( ) ~ t 的离散图2的观察,去掉时刻6及6.5以后的样本(这样作的原因见文章后面的评注),再利用最小二乘法进行拟合得lnpt( ) = 915158. -147577. t -147577. t其相关系数r = 0999887. ,由此得知pt() = 9429.33e 3 -147577. t我们作出pt() = 9429.33e 的图形,并将此图和图1放在一起图3,由图3及相关系数r = 0999887. 可以认为p()t 确实是负指数曲线-λtpt() ===Ae ,(A 9429.33,λ 147577. ).(2)及假设f，即h()0 = 0，解得kA --λtKtK -λ此式从数学上来看并不复杂，但要利用此式求出参数K和k却并非事，而参数K则需要在测试中使用，因此我们的问题归结为：如何利实际测量值和(2)及h(0)=0去决定参数k和K．这类数学问题称为参数识问题。

下面建立几个算法模型:算法模型Ⅰ.一般差分拟合法：将方程(2)离散化，记时间步长为T，利用前插公式得：hh-nn+1=-Kh +kpnnThKThk= ()1- + Tp ，(4)nn+1 n中hhtnTpptnT= (),()+ = + .nn002用差分法求解，其截断误差为oT()，显然大了些，为了提高精度准确度，最直接的方法是由插值方法获得更多的结点，缩短步长，使断误差减少。

如用三次样条插值法在每两个结点的中点进行插值，可2截断误差减少到原来的1/4，但仍然为oT()，且继续缩短步长，计算量将成倍增加。

算法模型Ⅱ.改进的差分拟合法：在这个算法中，我们注意方程(2)右端的线性项-Kh()t ，因此两边同Kt乘以e (积分因子)后可得: Ktde h()t Kt= kp()t e ，(5)dt对方程(5)利用差分离散化，并整理得：KTeh -hnn+1= kpnT即：--KT KThehkTep=+ ，(6)n+1 n n-KT 2此时截断误差为oe( T )，显然要比算法模型I 误差要小，同时若将(6)-KT -KT中的e 展开，即eKT=-1 +oT()，略去高阶无穷小，则得到：nn+1 n这恰好是方程(4)，由此可见利用积分因子后得到了一个比模型I 精度要高的一个算法模型。

对于离散方程(4)或(6)可以通过联立不同时刻的方程组求得一系列K值，但是由于在实际测量中存在随机误差，以及离散化的截断误差,使得这些K值不尽相同。

为了充分利用已测数据，我们利用最小二乘法拟合数据可得：hh= 0882488..+0078065p ，(7)nn+1 n在这里我们取t =1，步长T = 025. ，拟合的复相关系数r = 09999997.. 于0是将(7)与(4)式或(6)式比较可得参数K和k的值如下表所示：算法模型K的估计值k的估计值Ⅰ0.470048 0.31226Ⅱ0.5004 0.353841表1.算法Ⅰ算法Ⅱ的结果上述两个算法模型，计算简单，但对误差难以估计，并且对上述算法进行测试，两个算法对K具有稳定性，而对另一个参数k却不稳定，同时也看到算法Ⅱ优于算法Ⅰ，测试方法是预先假定一组K和k，按为未离散的公式(3)计算ht()在各时刻的值作为原始数据，再用差分公式和~ ~最小二乘法求出K和k ，将它们与原假定值作比较，测试的结果见表2：~ ~算法K k K kⅠ 1 2 0.8848 1.46852 1 1.3378 0.61699Ⅱ 1 2 0.999999 1.885622 1 1.93992 1.00071表2.测试情况使用求得的估计值K和k代入(3)式并作其连续图,然后与离散图作比较，同样可以看出模型Ⅱ优于模型Ⅰ(图4对应于模型Ⅰ,图5对应于模型Ⅱ)。

图4图5下面我们将给出另外一种算法对上述结果进行改进.算法摸型Ⅲ：线性迭代算法如果设已给K和k的预测值K 和k ,记0 0KK= +δ, 。

kk= +η0 0其中δ和η称为K和k相对于K 和k 的校正值(简称校正值)，将它们代0 0入(3)式并将右端关于δ和η展开成Taylor 级数，同时略去δ和η的二次及二次以上的项(即高阶无穷小项)，得到()kA0 +η --+λδtKt()0ht() = ()ee-K +-δλkA0 --λλtKt00ηA --tKt≈()()ee-+ ee-+K -λ K -λ00---Kt λt Kt00te ee-+δAk [ - ]2K0 -λ ()K0 -λ~= ht()利用理论值和实测数据误差的平方和最小的原则来选取δ和η,即选取δ和η使n~ 2Δ()hhtht=-[() ()]∑iii=1最小.利用最小二乘法求得δ和η后,较正K 和k 得0 0KK= +δ , kk= +η0 0将得到的新的参数K和k作为新的预测值，用同样的方法继续校正，直至δ和η足够小为止。

我们采用模型II的结果作为预测值,进行上述迭代程序得到的结论如表3所示:迭代次预测值K 预测值校正值δ校正值ηΔ()h数k初始值0.50004 0.3538411 0.50004 0.353841 0.004573 0.066023 66.1172 0.504613 0.419864 -0.000667 0.000025 65.3778-6 -63 0.503946 0.419889 -1302. ×10 -5459. ×10 65.3557-7 -74 0.503945 0.419884 -3053. ×10 -4199. ×10 65.3557结果值0.503945 0.419884表3.算法Ⅲ的迭代结果由表3可见算法模型Ⅲ的优越性与准确性，并且得到K和k的最佳拟合值为：K=0.503945, k=0.419884-7这种算法收敛速度很快，并且得到K值误差数量级为10四.模型的评价及注记(1)我们所建立的前两个算法模型计算简单,但是稳定性较差；第三个算法模型是稳定的，并且具有快速收敛性，可获得较精确的脑血流量系数K。