局部脑血流的测定 数学建模

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

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

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

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

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

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

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

若某受试者的测试数据如下:
时间(分) 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 ln2
1 2 3
= kp(),t = Kh(),t =- ht(),
dt dt dt τ
而
Δht() = Δh ()t -Δh ()t +Δh ()t ，
123
于是
dh dh dh dh
123
=-+ ,
dt dt dt dt
即
dh ln2
=-Kht()+kpt()- ht(). (
dt τ
133
由于在测试时放射性同位素(如Xe)的半衰期τ一般很大,而测试
时间又很短(大约十几分钟左右),由此假定τ→+∞,于是(1)式变为:
dh=-Kh()t +kp()t . (dt
2. 算法模型的建立与改进
在建立算法模型之前,首先必须对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 确实是负指数曲线-λt
pt() ===Ae ,(A 9429.33,λ 147577. ).
(2)及假设f，即h()0 = 0，解得
kA --λtKt
K -λ
此式从数学上来看并不复杂，但要利用此式求出参数K和k却并非
事，而参数K则需要在测试中使用，因此我们的问题归结为：如何利
实际测量值和(2)及h(0)=0去决定参数k和K．这类数学问题称为参数
识问题。

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

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

算法模型Ⅱ.改进的差分拟合法：
在这个算法中，我们注意方程(2)右端的线性项-Kh()t ，因此两边同Kt乘以e (积分因子)后可得: Ktde h()t Kt= kp()t e ，(5)
dt
对方程(5)利用差分离散化，并整理得：
KT
eh -h
nn+1
= kp
n
T
即：
--KT KT
hehkTep=+ ，(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.4685
2 1 1.3378 0.61699
Ⅱ 1 2 0.999999 1.88562
2 1 1.93992 1.00071
表2.测试情况
使用求得的估计值K和k代入(3)式并作其连续图,然后与离散图作
比较，同样可以看出模型Ⅱ优于模型Ⅰ(图4对应于模型Ⅰ,图5对应于模型Ⅱ)。

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

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

我们采用模型II的结果作为预测值,进行上述迭代程序得到的结论如表3所示:
迭代次预测值K 预测值校正值δ校正值ηΔ()h
数k
初始值0.50004 0.353841
1 0.50004 0.353841 0.004573 0.066023 66.117
2 0.50461
3 0.41986
4 -0.000667 0.00002
5 65.3778
-6 -6
3 0.503946 0.419889 -1302. ×10 -5459. ×10 65.3557
-7 -7
4 0.50394
5 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。

利用得到的结果，(3)式和已知的数据作ht()~ t 的连续图和离散图，如图6。

由图6显而易见我们的结果的精确度是非常高的。

图6
(2)在建模时忽略了同位数的衰变已及动脉血从肺部到脑部所需要
的时间，如在模型考虑这些因素后，只须在测试中测得这些因素的数值，用上述方法仍是容易实现的。

(3)在处理呼出气计数率的曲线时忽略了后面一些数据，我们认为这是合理的，因为任何测量仪器都有一定的精度要求，当呼出气中同位素计数率小到一定程度时，仪器是无法测出的，此时仪器的读数必将显示为零，故而读数是零并非说明计数率为零，所以不考虑后面为零的读数是合乎实际情况的。