1、局部脑血流的测定一. 问题简介 脑血流量是诊断和治疗脑梗塞,脑出血,动脉瘤和先天性动,静脉血管畸形等脑血管疾病的主要依据。测定脑血流量可为研究人脑在不同的病理和生理条件下的功能提供客观指标，它对研究脑循环药物的药理作用也很有帮助。所以人们长期致力于寻找有效地测定脑血流量的方法。 近年来出现了以放射性同位素作示踪剂测定人脑局部血流量的方法。这种方法大致可描述如下:由受试者吸入某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处放射性同位素的计数率(简称计数率) ，同时测量他呼出气的计数率。 由于动脉血将肺部的放射性同位素输送至大脑，使脑部同位素增加，而脑血流又将同位素带离，使
2、同位素减少，实验证明由脑血流引起局部地区计数率下降的速率与当时该处的记数率成正比，其比例系数反映了该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。动脉血从肺输送同位素至大脑引起脑部计数率上升的速率与当时呼出气的计数率成正比。 若某受试者的测试数据如下: 时 间(分) 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
3、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 呼出气记数
4、率 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、次仪器测量为相互独立事件，各测量值无记忆相关； 5. 放射性同位素在人体内传递是从吸入气体(含有放射物 )开始的，并假定一次吸入，因此认为同位素在肺中瞬时达到最大浓度； 6. 在吸入气体瞬时，脑中放射物记数率为零； 7. 脑血流量与脑血流系数 K 成单值函数关系，求得后者即可确定前者。三. 模型的建立与分析 由于已知脑局部同位素的增减与已测定的头部记数率 ht()和呼出气记数率 p()t 成正比关系,于是很自然地确定以脑部同位素量,即脑部记数率作为讨论对象.。 1.原始模型的建立 设某时刻 t 0 时,头部记数率为 ht(),在 t 时段后记数率 ht()+t ，由假定可知, 头部记数率的增量
6、 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由
7、于在测试时放射性同位素(如 Xe)的半衰期 一般很大,而测试时间又很短(大约十几分钟左右 ),由此假定 +,于是 (1)式变为: dh=-Kh()t +kp()t . (dt2. 算法模型的建立与改进 在建立算法模型之前,首先必须对 p()t 进行预测。作 p()tt 和-tln pt( ) t 的离散图(图 1 和图 2),由此发现 p()t 与 t 有近似于 Ae 的函数关系。 通过对 ln pt( ) t 的离散图 2 的观察,去掉时刻 6 及 6.5 以后的样本(这样作的原因见文章后面的评注),再利用最小二乘法进行拟合得 lnpt( ) = 915158. -147577. t -14
8、7577. 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 -tKtht() = (ee- ) ， (3) K -此式从数学上来看并不复杂，但要利用此式求出参数 K 和 k 却并非事，而参数 K 则需要在测试中使用，因此我们的问题归结为：如何利实际
9、测量值和(2)及 h(0)=0 去决定参数 k 和 K这类数学问题称为参数识问题。 下面建立几个算法模型: 算法模型.一般差分拟合法： 将方程(2)离散化，记时间步长为 T，利用前插公式得： hh-nn+1=-Kh +kp nnThKThk= ()1- + Tp ， (4) nn+1 n中 hhtnTpptnT= (),()+ = + . nn002用差分法求解，其截断误差为 oT()，显然大了些，为了提高精度准确度，最直接的方法是由插值方法获得更多的结点，缩短步长，使断误差减少。如用三次样条插值法在每两个结点的中点进行插值，可2 截断误差减少到原来的 1/4，但仍然为 oT()，且继续缩短步
10、长，计算量将成倍增加。 算法模型.改进的差分拟合法： 在这个算法中，我们注意方程(2)右端的线性项-Kh()t ，因此两边同 Kt 乘以 e (积分因子)后可得: Ktde h()t Kt= kp()t e ， (5) dt对方程(5)利用差分离散化，并整理得： KTeh -hnn+1= kp nT即： -KT KThehkTep=+ ， (6) n+1 n n-KT 2此时截断误差为 oe( T )，显然要比算法模型 I 误差要小，同时若将(6)-KT -KT中的 e 展开，即 eKT=-1 +oT()，略去高阶无穷小，则得到： hKThk= ()1- + Tp nn+1 n这恰好是方程(4
11、)，由此可见利用积分因子后得到了一个比模型 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
12、.5004 0.353841 表 1.算法算法的结果 上述两个算法模型，计算简单，但对误差难以估计，并且对上述算法进行测试，两个算法对 K 具有稳定性，而对另一个参数 k 却不稳定，同时也看到算法优于算法，测试方法是预先假定一组 K 和 k，按为未离散的公式(3)计算 ht()在各时刻的值作为原始数据，再用差分公式和 最小二乘法求出 K 和 k ，将它们与原假定值作比较，测试的结果见表 2： 算法 K k K k 1 2 0.8848 1.4685 2 1 1.3378 0.61699 1 2 0.999999 1.88562 2 1 1.93992 1.00071 表 2.测试情况 使用求得
13、的估计值 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 +-0kA0
14、 -tKt00A -tKt ()()ee-+ ee-+K - K -00-Kt t Kt00te ee-+Ak - 02K0 - ()K0 - = ht()利用理论值和实测数据误差的平方和最小的原则来选取 和 ,即选取 和 使 n 2()hhtht=-() () iii=1最小.利用最小二乘法求得 和 后,较正 K 和 k 得 0 0KK= + , kk= + 0 0将得到的新的参数 K 和 k 作为新的预测值，用同样的方法继续校正，直至 和 足够小为止。 我们采用模型 II 的结果作为预测值,进行上述迭代程序得到的结论如表 3 所示: 迭代次 预测值 K 预测值 校正值 校正值 ()h 数
15、k 初始值 0.50004 0.353841 1 0.50004 0.353841 0.004573 0.066023 66.117 2 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.50
16、3945, k=0.419884 -7这种算法收敛速度很快，并且得到 K 值误差数量级为 10 四.模型的评价及注记 (1)我们所建立的前两个算法模型计算简单,但是稳定性较差；第三个算法模型是稳定的，并且具有快速收敛性，可获得较精确的脑血流量系数 K。利用得到的结果，(3)式和已知的数据作 ht() t 的连续图和离散图，如图 6。由图 6 显而易见我们的结果的精确度是非常高的。图 6 (2)在建模时忽略了同位数的衰变已及动脉血从肺部到脑部所需要的时间，如在模型考虑这些因素后，只须在测试中测得这些因素的数值，用上述方法仍是容易实现的。 (3)在处理呼出气计数率的曲线时忽略了后面一些数据，我们认为这是合理的，因为任何测量仪器都有一定的精度要求，当呼出气中同位素计数率小到一定程度时，仪器是无法测出的，此时仪器的读数必将显示为零，故而读数是零并非说明计数率为零，所以不考虑后面为零的读数是合乎实际情况的。
