局部脑血流的测定一. 问题简介脑血流量是诊断和治疗脑梗塞,脑出血,动脉瘤和先天性动,静脉血管畸形等脑血管疾病的主要依据。
测定脑血流量可为研究人脑在不同的病理和生理条件下的功能提供客观指标，它对研究脑循环药物的药理作用也很有帮助。
所以人们长期致力于寻找有效地测定脑血流量的方法。
近年来出现了以放射性同位素作示踪剂测定人脑局部血流量的方法。
这种方法大致可描述如下:由受试者吸入某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处放射性同位素的计数率(简称计数率)，同时测量他呼出气的计数率。
由于动脉血将肺部的放射性同位素输送至大脑，使脑部同位素增加，而脑血流又将同位素带离，使同位素减少，实验证明由脑血流引起局部地区计数率下降的速率与当时该处的记数率成正比，其比例系数反映了该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。
动脉血从肺输送同位素至大脑引起脑部计数率上升的速率与当时呼出气的计数率成正比。
试建立确定脑血流系数的数学模型并计算上述受试者的脑血流系数。
备注:该题目是市(1990 年)大学生数学建模竞赛A 题。
二.模型的假定= − 1. 脑部计数率(记为 h (t ) )的上升只与肺部的放射性同位素有关，上 升速度与呼出气的记数率(记为 p (t ) )成正比，比例系数记为 k ；2. 脑部记数率 h (t ) 的下降只与该处脑血流量有关，其下降速度正比 于 h (t ) ，比例系数为脑血流系数，记为 K ，这里忽略了放射性元素的衰 变和其它因素；3. 脑血流量在测定期间恒定，心脏博动，被测试者大脑活动，情感 波动等带来的变化忽略不予考虑；4. 每次仪器测量为相互独立事件，各测量值无记忆相关；5. 放射性同位素在人体传递是从吸入气体(含有放射物)开始的， 并假定一次吸入，因此认为同位素在肺中瞬时达到最大浓度；6. 在吸入气体瞬时，脑中放射物记数率为零；7. 脑血流量与脑血流系数 K 成单值函数关系，求得后者即可确定前 者。
三. 模型的建立与分析 由于已知脑局部同位素的增减与已测定的头部记数率 h (t ) 和呼出气 记数率 p (t ) 成正比关系,于是很自然地确定以脑部同位素量,即脑部记数 率作为讨论对象.。
1.原始模型的建立 设某时刻 t ≥ 0 时,头部记数率为 h (t ) ,在 Δt 时段后记数率 h (t + Δt ) ， 由假定可知, 头部记数率的增量 Δh = h (t + Δt ) − h (t ) 仅与三个因素有关： (i) 肺动脉血将肺部的放射性同位素送至大脑,使脑部记数率增量 为 Δh 1 ；(ii) 脑血流将同位素带离,脑记数率下降为 Δh 2 ；(iii)放射性同位素自身有衰减引起记数率下降量为 Δh 3 ,设其半衰期 为τ .因此,由医学试验及假定有dh 1 = kp (t ), dh 2= Kh (t ),dh 3 ln 2h (t ),τdtdtdt而Δh (t ) = Δh 1 (t ) − Δh 2 (t ) + Δh 3 (t ) ，于是dh = dh 1 − dh 2 + dh 3 , dt dt dt dt即dh = − K h (t ) + kp (t ) − ln 2 h (t ) . (1)dt τ由于在测试时放射性同位素(如133 Xe )的半衰期τ 一般很大,而测试时间又很短(大约十几分钟左右),由此假定τ → +∞ ,于是(1)式变为:dh= − K h (t ) + kp (t ) . (2) dt2. 算法模型的建立与改进在建立算法模型之前,首先必须对 p (t ) 进行预测。
作 p (t ) ~ t 和ln p (t ) ~ t 的离散图(图 1 和图 2),由此发现 p (t ) 与 t 有近似于 Ae − λt 的函 数关系。
通过对 ln p (t ) ~ t 的离散图 2 的观察,去掉时刻 6 及 6.5 以后的样本点 (这样作的原因见文章后面的评注),再利用最小二乘法进行拟合得图1 图2ln p (t ) = 9.15158 − 1.47577t其相关系数 r = 0.999887 ,由此得知 p (t ) = 9429.33e −1.47577t .图3我们作出p(t ) = 9429.33e −1.47577t 的图形,并将此图和图1 放在一起见图3,由图3 及相关系数r = 0.999887 可以认为p(t )确实是负指数曲线p(t ) = Ae −λt , (A = 9429.33, λ= 1.47577) .由(2)及假设f，即h(0) = 0 ，解得h(t ) =kAK −λ(e −λt −e −Kt )，(3)此式从数学上来看并不复杂，但要利用此式求出参数K 和k 却并非易事，而参数K 则需要在测试中使用，因此我们的问题归结为：如何利用实际测量值和(2)及h(0)=0 去决定参数k 和K ．这类数学问题称为参数辨识问题。
下面建立几个算法模型:算法模型Ⅰ.一般差分拟合法：将方程(2)离散化，记时间步长为T，利用前插公式得：h −hn+1 n = −K hn + kp nT即h n+1 = (1 −KT )h n + kTp n ，(4)其中h n = h(t0 + nT ), p n = p(t0 + nT ).用差分法求解，其截断误差为o(T 2 ) ，显然大了些，为了提高精度和准确度，最直接的方法是由插值方法获得更多的结点，缩短步长，使截断误差减少。
如用三次样条插值法在每两个结点的中点进行插值，可使截断误差减少到原来的1/4，但仍然为o(T 2 ) ，且继续缩短步长，计算量将成倍增加。
算法模型Ⅱ.改进的差分拟合法：在这个算法中，我们注意方程(2)右端的线性项−Kh(t ) ，因此两边同乘以e Kt (积分因子)后可得:Ktde h(t)=kp(t )e Kt，(5)dt对方程(5)利用差分离散化，并整理得：e KT h n+1 −h nT= kp n即：h n+1 = e −KT h n + kTe −KT p n ，(6)此时截断误差为o(e −KT T 2 ) ，显然要比算法模型I 误差要小，同时若将(6)中的e −KT 展开，即e −KT= 1 −KT + o(T )，略去高阶无穷小，则得到：h n+1 = (1 −KT )h n + kTp n这恰好是方程(4)，由此可见利用积分因子后得到了一个比模型I 精度要高的一个算法模型。
对于离散方程(4)或(6)可以通过联立不同时刻的方程组求得一系列K 值，但是由于在实际测量中存在随机误差，以及离散化的截断误差, 使得这些K 值不尽相同。
为了充分利用已测数据，我们利用最小二乘法拟合数据可得：h n+1 = 0.882488h n + 0.078065 p n ，(7)在这里我们取t= 1 ，步长T = 0.25 ，拟合的复相关系数r = 0.9999997 . 于是将(7)与表1.算法Ⅰ算法Ⅱ的结果上述两个算法模型，计算简单，但对误差难以估计，并且对上述算法进行测试，两个算法对K 具有稳定性，而对另一个参数k 却不稳定，同时也看到算法Ⅱ优于算法Ⅰ，测试方法是预先假定一组K 和k ，按为未离散的公式(3)计算h(t ) 在各时刻的值作为原始数据，再用差分公式和最小二乘法求出~和~，将它们与原假定值作比较，测试的结果见表2：K kh表 2.测试情况使用求得的估计值 K 和 k 代入(3)式并作其连续图,然后与离散图作比较，同样可以看出模型Ⅱ优于模型Ⅰ(图 4 对应于模型Ⅰ,图 5 对应于模型Ⅱ)。
图 4图 5下面我们将给出另外一种算法对上述结果进行改进. 算法摸型Ⅲ：线性迭代算法如果设已给 K 和 k 的预测值 K 0 和 k 0 ,记K = K 0 + δ , 。
k = k 0 + η其中δ 和η 称为 K 和 k 相对于 K 0 和 k 0 的校正值(简称校正值)，将它们代 入(3)式并将右端关于δ 和η 展开成 Taylor 级数，同时略去δ 和η 的二次 及二次以上的项(即高阶无穷小项)，得到h (t ) = ( k 0 + η) A (e − λt − e −( K 0 +δ )t )K 0 + δ − λ≈ k 0 A K 0 − λ(e − λt − e − K 0t ) +ηAK 0 − λ(e − λt − e − K 0t ) + − K t− λt− K t+ δAk [ te 0 − e − e 0 ]0 K − λ ( K − λ ) 2 0 0 = ~(t )利用理论值和实测数据误差的平方和最小的原则来选取δ 和η ,即选取δ 和η 使 Δ(h ) = n ~ [h (t ) − (t )] ∑ ii =1h i 2 最小.利用最小二乘法求得δ 和η 后,较正 K 0 和 k 0 得K = K 0 + δ , k = k 0 + η将得到的新的参数 K 和 k 作为新的预测值，用同样的方法继续校正，直 至δ 和η 足够小为止。
我们采用模型 II 的结果作为预测值,进行上述迭代程序得到的结论如表 3.算法Ⅲ的迭代结果由表 3 可见算法模型Ⅲ的优越性与准确性，并且得到 K 和 k 的最佳 拟合值为：K =0.503945, k =0.419884这种算法收敛速度很快，并且得到 K 值误差数量级为10−7 . 四.模型的评价及注记(1)我们所建立的前两个算法模型计算简单,但是稳定性较差；第三 个算法模型是稳定的，并且具有快速收敛性，可获得较精确的脑血流量 系数 K 。
利用得到的结果，(3)式和已知的数据作 h (t ) ~ t 的连续图和离散 图，如图 6。
由图 6 显而易见我们的结果的精确度是非常高的。
图6 (2)在建模时忽略了同位数的衰变已及动脉血从肺部到脑部所需要的时间，如在模型考虑这些因素后，只须在测试中测得这些因素的数值，用上述方法仍是容易实现的。
(3)在处理呼出气计数率的曲线时忽略了后面一些数据，我们认为这是合理的，因为任何测量仪器都有一定的精度要求，当呼出气中同位素计数率小到一定程度时，仪器是无法测出的，此时仪器的读数必将显示为零，故而读数是零并非说明计数率为零，所以不考虑后面为零的读数是合乎实际情况的。