脑部血流测定 (2)

p(t)拟 合 前 后 对 比
2000
1500
1000
500
0
1
2
3
4
5
6
7
8
9
10
图2
我们作出 p(t) 9920 .9e1.4953t 的图形，并将其与原始数据放在一起对比，由
图 2 可以认为 p(t) 确实是负指数曲线
p(t) Aet , ( A 9920 .9, 1.4953 )
由（2）及假设知，即 h(0) 0 ，解得
5.1 模型一：一般差分拟合法 将方程（2）离散化，记时间步长为T ，利用前差分公式得： hn 1 hn Khn kpn T
即
其中
hn 1 (1 KT )hn kTpn
hn h(t0 nT ) pn p(t0 nT )
（4）
用差分法求解，其截断误差为 o(T 2 ) ，显然大了一些，为了提高精度和准确 度，最直接的方法是由插值获得更多的结点，缩短步长，使截断误差减少。如用 三次样条插值法在每两个结点的中点进行插值，可使截断误差减少到原来的 1 ，
2.5
1500
2
1000
1.5
1
500 0.5
0
0
2
4
6
8
10
0
0
2
4
6
8
10
图1
通过对 Ln( p(t)) ~ t 的离散图的观察，去掉时刻 6 及以后的样本点，再利用最
小二乘法进行拟合得
Ln( p(t)) 1.4953t 9.2024
由此可知
p(t) 9920 .9e1.4953t
2500
4 但仍然为 o(T 2 ) ，且继续缩短步长，计算量将成倍增加。 5.2 模型二：改进的差分拟合法
在这个算法中，我们注意方程（2）右端的线性项 Kh(t) ，因此两边同时乘
以 eKt （积分因子）后可得：
deKt h(t) kp(t)eKt dt
对方程（5）利用差分离散化，并整理得：
即：
eKT hn 1 hn kpn T
1.2 问题
用放射性同位素测定大脑局部血流量的方法如下：有受试者吸入含有某种放 射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处的放 射性记数率（简称记数率），同时测量他呼出气的计数率。
由于动脉血将肺部的放射性同位素传送至大脑，使脑部同位素增强，而脑血 流又将同位素带离，是同位素减少。实验证明由脑血流引起局部地区记数率下降 的速率与当时该处的记数率成正比。其比例系数反应该处的脑血流量，被称为脑 血流量系数，只要确定该系数即可推算出脑血流量。动脉血从肺输送同位素至大 脑引起脑部记数率上升的速率与当时呼出气的记数率成正比。若某受试者的测试 数据如附录一所示 根据以上题目所给的条件及数据，回答以下问题： 问题一：建立确定脑部血流系数的数学模型； 问题二：计算上述受试者的脑血流系数。
（即高阶无穷小项），得到
h(t) (k 0 ) A (et e(K 0 )t ) K0
k 0 A (et eK 0t ) A (et eK 0t )
K0
K0
Ak0[ teK 0t K0
et eK 0t (K0 )2
]
h~(t)
利用实际理论值和实测数据误差的平方和最小的原则来选取 和 ，即选取
（1） 肺动脉血将肺部的放射性同位素送至大脑，使脑部记数率增量 h1 ； （2） 脑血流将同位素带离，脑记数率下降减量 h2 ； （3） 放射性同位素自身有衰减引起记数率下降剑量为 h3 ，设其半衰期为
。
因此，由医学试验及假定有
dh1 kp(t) dt
dh2 Kh(t) dt
dh3 ln 2 h(t) dt
hn 1 eKT hn kTeKT pn
（5） （6）
此时截断误差为 o(eKTT 2 ) ，显然要比算法模型一误差小，同时若将（6）中
的 eKT 展开，即 eKT 1 KT o(T ) ，略去高阶无穷小，则得到：
hn 1 (1 KT )hn kTpn
这恰好是方程（4），由此可见利用积分因子后得到了一个比模型一精度要高的 一个算法模型。
2. 脑部计数率的上升只是由于肺将同位素通过血液循环输送至大脑，其上升速 率与当时呼出气体的计数率成正比。
3. 在实验开始瞬间，脑中无放射物，此时头部计数率为零。 4. 假设脑血流恒定，实验中不随人的意志的变化而变化。
5. 由与实验时常用的放射性同位素的半衰期很大，所以不考虑同位素衰变对实 验的影响。
四、 符号说明
ht
表示 t 时刻的脑部计数率
pt
呼出气的记数率
k
脑血流量系数
脑部记数率上升的速率与当时呼出气体记
K
数率之比
肺动脉血将肺部的放射性同位素送至大脑，
h1
使脑部记数率增量
h2
脑血流将同位素带离，脑记数率下降减量
放射性同位素自身有衰减引起记数率 模型建立及求解
和 使
(h) n [h(ti) h~(ti)]2 i 1
最小。利用最小二乘法求得 和 后，校正 K0 和 k 0 得
K K0 k k0
将得到的新的参数 K 和 k 作为新的预测值，用同样的方法继续校正，直到 和
足够小为止。 我们采用模型二的结果作为预测值，进行上述迭代程序得到的结论如下表所
本文在科学分析题目所提供的相关数据的基础上，首先科学预测并拟合得出 呼出气记数率的变化近似于 9920.9e1.4953 ，然后由浅至深建立了以最先而成拟合 法为主的三个数学模型，三个模型联系紧密，都是在上一个模型的基础上精确化 得到的。
最后得到 ht 4201 .9531 e e 1.4953 0.503945 可获得较精确的脑部记数率上升
而
h(t) h1(t) h2(t) h3(t)
于是 即
dh dh1 dh2 dh3 dt dt dt dt
dh kp(t) Kh(t) ln 2 h(t)
dt
（1）
由于在测试时放射性同位素（如 133 Xe ）的半衰期 一般很大，而测试时间
又很短（大约十几分钟），由此假定 ，于是式（1）变为：
示：
迭代次数 预测值 K 预测值 k
校正值
校正值
(h)
初始值 0.5000409 0.3538404
1
0.5000409 0.3538404 0.0800260 0.0009089 5.904148
2
0.5800669 0.3547493 -0.000667 0.000025 6.53778
3
0.503946 0.419889 1.302106 5.459106 6.53557
对于离散方程（4）或（6）可以通过联立不同时刻的方程组求得一系列 K 值， 但是由于在实际测量中存在随机误差，以及离散化的截断误差，使得这些 K 值不 尽相同。为了充分利用已测数据，我们利用最小二乘发拟合数据可得
hn 1 0.8825hn 0.0781pn
(7)
在这里我们取 t0 1，步长T 0.25，拟合得负相关系数 r 0.9999997，最大绝对 误差为 0.6931544。于是（7）与（4）式或（6）式比较可得参数 K 和 k 的值如下 表所示：
我们参赛选择的题号是（从 A/B/C/D 中选择一项填写）：
我们的参赛报名号为（如果赛区设置报名号的话）：
所属学校（请填写完整的全名）：
西安理工大学
参赛队员 (打印并签名) ：1.
何吉
2.
余蓉
3.
张苗
指导教师或指导教师组负责人 (打印并签名)：
邹学文
日期： 2011 年 8 月 10 日
2011 高教社杯全国大学生数学建模竞赛
全国大学生数学建模竞赛
承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则. 我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨 询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。 我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料 （包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中 明确列出。 我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的 行为，我们将受到严肃处理。
h(t) kA (et eKt ), K
（3）
此式从数学上来看并不复杂，但要利用此式求出参数 K 和 k 却并非易事，而
参数 K 则需要在测试中使用，因此我们的问题归结为：如何利用实际测量值和（2）
及 h(0) 0 去决定参数 K 和 k 。这类数学问题称为参数辨识问题。
三、 模型假设
1. 脑部计数率的下降只是由于脑血流将同位素带离脑部，其下降速率与头部计 数率成正比。
在各时刻的值作为原始数据，再用差分公式和最小二乘法求出 K~ 和 k~ ，将它们与
原假定值作比较，测试结果如下表 3：
表 3 测试结果
K
k
K~
k~
模型一
2
1
1.5675099
0.5930288
1
2
0.8833142
1.3511987
模型二
2
1
1.9895155
0.9751798
1
2
0.9980967
400
200
200
0
0
0
5
10
0
5
10
图 3 两模型对比图
5.3 模型三：线性迭代算法
如果设已给出 K 和 k 的预测值 K0 和 k 0 ，记 K K0 k k0
其中 和 称为 K 和 k 相对于 K0 和 k 0 的校正值（简称校正值），将它们代入（3）