脑部血流测定-(2)

全国大学生数学建模竞赛承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）：我们的参赛报名号为（如果赛区设置报名号的话）：所属学校（请填写完整的全名）：西安理工大学参赛队员(打印并签名)：1.何吉2.余蓉3.张苗指导教师或指导教师组负责人(打印并签名)：邹学文日期：2011年8月10日2011高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：局部脑血流测定摘要脑血流量是诊断和治疗脑梗塞、脑出血、动脉瘤和先天性动静脉血管畸形等脑血管疾病的主要依据。测定脑血流量可以为研究人脑在不同的病理和生理条件下的功能提供客观指标，它对研究脑循环药物的药理作用有很大帮助。本题要求建立合理的脑血流系数模型及求解题目给数据的脑血流系数。本文在科学分析题目所提供的相关数据的基础上，首先科学预测并拟合得出呼出气记数率的变化近似于4953.19.9920e，然后由浅至深建立了以最先而成拟合法为主的三个数学模型，三个模型联系紧密，都是在上一个模型的基础上精确化得到的。最后得到503945.04953.19531.4201eeth可获得较精确的脑部记数率上升的速率与当时呼出气体记数率之比值与脑血流量系数503945.0K和419884.0k。关键词：脑血流测定最小二乘法拟合差分线性迭代一、问题重述1.1问题背景在发达国家中，心脑血管疾病是威胁人们生命的最主要疾病之一。在我国，由于人民生活的改善和健康水准的提高，其他疾病的发病率下降，防治水平提高，心脑血管的发病率及其导致的死亡率却相对地上升了。脑血流量是诊断和治疗脑梗塞，脑出血，动脉瘤和先天性动脉和静脉血管畸形等脑血管疾病的主要依据。测量脑血流量可为研究人脑在不同的病理和生物条件下（入脑外伤，脑循环停顿，缺氧等）的功能提供客观指标，它对研究脑循环药物的药理作用也很有帮助。所以人们长期致力于寻找有效地测定脑血流的方法。1.2问题用放射性同位素测定大脑局部血流量的方法如下：有受试者吸入含有某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处的放射性记数率（简称记数率），同时测量他呼出气的计数率。由于动脉血将肺部的放射性同位素传送至大脑，使脑部同位素增强，而脑血流又将同位素带离，是同位素减少。实验证明由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比。其比例系数反应该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比。若某受试者的测试数据如附录一所示根据以上题目所给的条件及数据，回答以下问题：问题一：建立确定脑部血流系数的数学模型；问题二：计算上述受试者的脑血流系数。二、问题分析由于已知脑部局部同位素的增减与已测定的头部记数率)(th和呼出气记数率)(tp成正比关系，于是很自然地确定以脑部同位素量，即脑部记数率作为讨论对象。2.1原始问题的分析设某时刻0t，头部记数率为)(th，在t时段后记数率为)(tth，由假定可知，头部记数率的增量)()(thtthh仅与三个因素有关：（1）肺动脉血将肺部的放射性同位素送至大脑，使脑部记数率增量1h；（2）脑血流将同位素带离，脑记数率下降减量2h；（3）放射性同位素自身有衰减引起记数率下降剑量为3h，设其半衰期为。因此，由医学试验及假定有)(2ln)()(321thdtdhtKhdtdhtkpdtdh而)()()()(321thththth于是dtdhdtdhdtdhdtdh321即)(2ln)()(thtKhtkpdtdh（1）由于在测试时放射性同位素（如Xe133）的半衰期一般很大，而测试时间又很短（大约十几分钟），由此假定，于是式（1）变为：)()(tKhtkpdtdh（2）2.2算法的建立和改进在建立算法模型之前，首先必须对)(tp进行预测。作ttp~)(和ttpLn~))((的离散图（图1），由此发现)(tp与t有近似于tAe的函数关系。024681005001000150020002500p(t)024681000.511.522.533.5Ln(p(t))图1通过对ttpLn~))((的离散图的观察，去掉时刻6及以后的样本点，再利用最小二乘法进行拟合得2024.94953.1))((ttpLn由此可知tetp4953.19.9920)(1234567891005001000150020002500p(t)拟合前后对比图2我们作出tetp4953.19.9920)(的图形，并将其与原始数据放在一起对比，由图2可以认为)(tp确实是负指数曲线)4953.1,9.9920(,)(AAetpt由（2）及假设知，即0)0(h，解得),()(KtteeKkAth（3）此式从数学上来看并不复杂，但要利用此式求出参数K和k却并非易事，而参数K则需要在测试中使用，因此我们的问题归结为：如何利用实际测量值和（2）及0)0(h去决定参数K和k。这类数学问题称为参数辨识问题。三、模型假设1.脑部计数率的下降只是由于脑血流将同位素带离脑部，其下降速率与头部计数率成正比。2.脑部计数率的上升只是由于肺将同位素通过血液循环输送至大脑，其上升速率与当时呼出气体的计数率成正比。3.在实验开始瞬间，脑中无放射物，此时头部计数率为零。4.假设脑血流恒定，实验中不随人的意志的变化而变化。5.由与实验时常用的放射性同位素的半衰期很大，所以不考虑同位素衰变对实验的影响。四、符号说明th表示t时刻的脑部计数率tp呼出气的记数率k脑血流量系数K脑部记数率上升的速率与当时呼出气体记数率之比1h肺动脉血将肺部的放射性同位素送至大脑，使脑部记数率增量2h脑血流将同位素带离，脑记数率下降减量3h放射性同位素自身有衰减引起记数率下降减量放射性元素半衰期五、模型建立及求解5.1模型一：一般差分拟合法将方程（2）离散化，记时间步长为T，利用前差分公式得：nnnnkpKhThh1即nnnkTphKTh)1(1（4）其中)()(00nTtppnTthhnn用差分法求解，其截断误差为)(2To，显然大了一些，为了提高精度和准确度，最直接的方法是由插值获得更多的结点，缩短步长，使截断误差减少。如用三次样条插值法在每两个结点的中点进行插值，可使截断误差减少到原来的41，但仍然为)(2To，且继续缩短步长，计算量将成倍增加。5.2模型二：改进的差分拟合法在这个算法中，我们注意方程（2）右端的线性项)(tKh，因此两边同时乘以Kte（积分因子）后可得：KtKtetkpdtthde)()(（5）对方程（5）利用差分离散化，并整理得：nnnKTkpThhe1即：nKTnKTnpkTeheh1（6）此时截断误差为)(2TeoKT，显然要比算法模型一误差小，同时若将（6）中的KTe展开，即)(1ToKTeKT，略去高阶无穷小，则得到：nnnkTphKTh)1(1这恰好是方程（4），由此可见利用积分因子后得到了一个比模型一精度要高的一个算法模型。对于离散方程（4）或（6）可以通过联立不同时刻的方程组求得一系列K值，但是由于在实际测量中存在随机误差，以及离散化的截断误差，使得这些K值不尽相同。为了充分利用已测数据，我们利用最小二乘发拟合数据可得nnnphh0781.08825.01(7)在这里我们取10t，步长25.0T，拟合得负相关系数9999997.0r，最大绝对误差为0.6931544。于是（7）与（4）式或（6）式比较可得参数K和k的值如下表所示：表2模型一和二结果Kk模型一0.47004850.3122599模型二0.50004090.3538404上述两种算法模型，计算简单，但对误差难以估计，并且对上述算法进行测试，两个算法对K具有稳定性，而另一种参数k却不稳定，同时也看到算法二优越于算法一，测试方法是预先假定一组K和k，按为未离散的公式（3）计算)(th在各时刻的值作为原始数据，再用差分公式和最小二乘法求出K~和k~，将它们与原假定值作比较，测试结果如下表3：表3测试结果KkK~k~模型一21121.56750990.88331420.59302881.3511987模型二21121.98951550.99809670.97517981.7341481使用求得的估计值K和k代入（3）式并作出其连续图，然后与离散图作比较，同样可以看出模型二优于模型一。（如图3）051002004006008001000120014001600模型一原始数据算得结果对应数据051002004006008001000120014001600模型二原始数据算得结果对应数据图3两模型对比图5.3模型三：线性迭代算法如果设已给出K和k的预测值0K和0k，记00kkKK其中和称为K和k相对于0K和0k的校正值（简称校正值），将它们代入（3）式并将右端关于和展开成泰勒级数，同时略去和的二次及二次以上的项（即高阶无穷小项），得到)(~])([)()()()()(2000000)(0000000thKeeKteAkeeKAeeKAkeeKAkthtKttKtKttKttKt利用实际理论值和实测数据误差的平方和最小的原则来选取和，即选取和使niiiththh12)](~)([)(最小。利用最小二乘法求得和后，校正0K和0k得00kkKK将得到的新的参数K和k作为新的预测值，用同样的方法继续校正，直到和足够小为止。我们采用模型二的结果作为预测值，进行上述迭代程序得到的结论如下表所示：迭代次数预测值K预测值k校正值校正值)(h初始值0.50004090.353840410.50004090.35384040.08002600.00090895.90414820.58006690.3547493-0.0006670.0000256.5377830.5039460.419889610302.1610459.56.5355740.5039450.419884710053.3710199.46.53557结果值0.5039450.419884由表3可见算法模型三的优越性与准确性，并且得到K和k的最佳拟合值为：419884.0503945.0kK这种算法收敛速度很快，并且得到K值误差数量级为710。最终将得到的K，k代入(3)式),()(KtteeKkAth得到脑部计数率的数学模型503945.04953.19531.4201eeth模型三的拟合结果与原数据比较如图41234567891002004006008001000120014001600模型三原始数据算得结果对应数据图4模型三拟合图六、模型评价与推广6.1模型的评价我们所建立的前两个算法模型计算简单，但是稳定性较差；第三个算法模型较为稳定，并且具有快速收敛性，可获得较精确的脑血流量系数K。在处理呼出气记数率的曲线时忽略了后面一些数据，我们认为这是合理的，因为任何