局部脑血流测定欧阳家百（2021.03.07）摘要本文主要对人体大脑局部脑血流量进行测定，实验使受试者吸入某种放射性同位素的气体，定时测量放射性计数率和呼出气的计数率，由计数率变化速率与计数率和呼出气计数率的关系，求解头部计数率的随时间变化的关系。
针对问题1，首先根据题设可知：由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比与动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比的两个关系，得到脑部计数率的变化量的二元一阶线性非齐次常微分方程：；采用消元法，引入呼出气记数率与时间的关系函数，设定初始值：，可建立一阶线性非齐次常微分方程模型：，进行求解。
针对问题2，对上述模型进行求解，首先对原始数据脑部计数率与时间，呼出气计数率与时间的关系用进行拟合，得到拟合曲线，由曲线看出呼出气计数率与时间大致成指数关系，进而对呼出气计数率进行取对数的数据变化，用进行一次多项式拟合，拟合结果得到：。
将带入微分方程根据一阶线性非齐次常微分方程的通解得。
用MATLAB对其进行最小二乘法拟合，求得正比系数，。
问题二结果检验：1、初值检验：将带入，得与所给初始值1534近似相等，误差非常小，验证了结果的准确性；2、差值检验：由图得差值在直线上下波动较小。
因此结果比较准确。
关键字脑血流量系数常微分方程模型最小二乘法差值图一．问题重述用放射性同位素测定大脑局部血流量的方法如下：由受试者吸入含有某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处的放射性记数率（简称记数率），同时测量他呼出气的记数率。
由于动脉血将肺部的放射性同位素传送至大脑，使脑部同位素增加，而脑血流又将同位素带离，使同位素减少。
实验证明由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比。
其比例系数反应该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。
动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比。
某受试者的测试数据见附表1。
根据题目所给条件与数据，求解一下问题：1.建立确定脑部血流系数的数学模型；2.计算上述受试者的脑血流系数。
二．问题分析2.1 对问题1的分析：针对问题1，题目中给出了动脉血，脑血流对脑部计数率的影响。
首先，脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比，且比例系数反应该处的脑血流量。
另外，脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比。
由这两个正比关系即可得到脑部地区计数率总的变化率与时间的关系，列出微分方程，建立微分方程数学模型。
2.2 对问题2的分析：针对问题2，由问题1建立的微分方程模型进行求解。
考虑模型是二元一阶方程，无法求解。
我们对呼出气的计数率与时间的数据进行处理，用matlab进行拟合得到它们之间的关系方程，带入模型，模型变为一阶线性常微分方程，进而可以求解。
三．模型假设1.假设题目所给数据均真实可靠；2.假设受试者的脑血流量不受吸入放射性同位素气体的影响；3.假设受试者在吸入放射性同位素气体前，脑中无这种放射性同位素气体；4.假设脑部计数率的下降只与脑血流有关，且下降速率与该处的计数率成正比；5.假设脑部计数率的上升只与动脉血有关，且上升速率与当时呼出气的计数率成正比；6.假设每次测量的数据均是相互独立的。
四．符号说明符号意义表示时间时刻头部计数率时刻呼出气计数率脑部计数率下降的速率与该处计数率成正比关系的比例系数脑部计数率上升的速率与当时呼出气的计数率成正比关系的比例系数自定义常数，误差的大小差值五．模型的建立与求解5.1.1建模准备过程分析：以脑部计数率为研究对象，脑部计数率的变化分两个过程：1、脑血流使得脑部计数率下降，并且下降速率与该时刻脑部计数率成正比；2、动脉血使得头部计数率上升，并且上升速率与该时刻呼出气计数率成正比。
如图1：图1 头部计数率变化流程图 5.1.2 建模过程根据头部计数率变化流程图建立以下模型：设时刻头部计数率为，呼出气计数率为，经过时刻，由脑血流引起的头部计数率的变化，；由动脉血引起的头部计数率的变化，则经过时刻头部计数率的总变化量，即：，此方程为二元一阶常系数线性常微分方程。
动脉血头部计数率 脑血流上升：速率与该时刻头部计数率成正比下降：速率与该时刻呼出气计数率成正比消元法求解：该方程为二元方程，不能求解，考虑消去。
引入呼出气计数率与时间的函数关系：，带入原方程得：，即：此方程为一阶线性非齐次常微分方程。
设定初始值：，即求解：5.2.1 模型求解受试者脑血流系数的计算：将原始数据脑部计数率与时间，呼出气计数率与时间的关系用matlab 进行拟合，得到拟合曲线如图2，图2 计数率随时间变化趋势图由图可以看出呼出气计数率与时间大致呈指数函数关系，因此，对呼出气数据进行取对数变换，得表2： 表2 呼出气计数率对数变换表取对数大于0的部分，用MATLAB 进行一次多项式拟合，得拟合系数，拟合曲线如图3：取对数 7.7102 7.3356 6.9603 6.5848 6.2106 5.8348 5.4596 5.0876 4.7095 4.3307 3.9512 3.5835 3.2189 2.8332 2.4849 2.0794 1.7918 1.3863 1.0986 0.6931图3 对数变换一次拟合直线与原始数据得到很好的匹配。
取对数后，即，5.2.2 残差分析：残差平方和的概念：为了明确解释变量和随机误差各产生的效应是多少，统计学上把公式数据点与它在回归直线上相应位置的差异称残差，把每个残差平方后加起来称为残差平方和。
对所求的的函数进行数据残差分析：用MATLAB工具求得该残差平方和为：，残差平方和很小，说明误差很小。
5.2.3数据检验：绘制原始数据与函数的对比图，如图4：图4 数据检验图将代入原方程：根据线性一阶非齐次微分方程的通式及其通解形式，解得：令得，其中。
5.2.4模型结果：采用最小二乘法进行拟合，拟合曲线见图5：图5 最小二乘拟合曲线得到参数：,根据求得动脉血头部计数率上升系数5.2.5模型检验：5.2.5.1初值分析检验当时，代入得：时，与所给初始值1534近似相等所得误差为：误差非常的小，因此验证了该模型的准确性。
5.2.5.2差值分析检验设时刻头部计数率的真实值为表示，拟合值为，差值。
做时间——差值图：图6 差值分析图由图可以看出，差值在直线上下波动，起伏很小，验证了结果的准确性。
六．模型评价与推广6.1.模型的评价模型的优点：模型属于微分方程模型，比较简单，但结果比较准确。
模型多次利用MATLAB进行数据拟合，且拟合结果均与实际相符合，对呼出气计数率与时间的关系先进行拟合，再提出猜想，最后进行验证，证明正确性。
模型求解采用最小二乘法拟合，最后将结果做差值图进行验证，得出较小的误差与分析，由此可以看出，模型结果比较准确，与实际相符合。
因此模型对实际脑部血流量的测定有很好的指导意义。
模型的缺点：本模型在建立的过程中没有考虑这种放射性同位素的衰变，以及动脉血从肺部到脑部所需要的时间，因此结果比较理想化，可能与实际存在一定误差。
6.2.模型的推广本模型可以推广到其他用放射性同位素测试的实际问题中，找出所研究问题与可以放射性同位素之间的关系，同样列出常微分方程模型进行求解。
同时该模型在医疗方面，可对病人病情进行检测。
具有很好的实际指导意义。
当考虑同位素的衰变，动脉血从肺部到脑部所需要的时间等因素后，可以实际测得这些数据，用本模型依然可以实现。
七．参考文献[1]曹卫华,郭止.最优化技术方法及MATLAB的实现[M],北京：化学工业出版社，2005.1[2]王家文，王皓，刘海.MATLAB7.0编程基础[M],北京：机械工业出版社，2005.7[3]刘志平，石林英.最小二乘法原理及其MATLAB实现[J],中国科技西部，2008，17(7):33-34八．附录表1某受试者的测试数据附源程序代码：%data.mclear;close all;clc;a=xlsread('C:\Users\谷柏辰\Desktop\data.xls')t=1:0.25:10;plot(t,a(1,:),'-b');hold onplot(t,a(2,:),'-r');function f=fun(x,xdata)n=length(xdata);for i=1:nf(i)=x(1)*exp(-1.4808*xdata(i))+x(2)*exp(-x(3)*xdata(i));endclcclear;t=[1:0.25:5.75];m=[2231 1534 1054 724 498 342 235 162 111 76 52 36 25 17 12 8 6 4 3 2];plot(t,m,'*');hold on;y=exp(-1.4808*t+9.1648);plot(t,y,'r');title('M函数原始数据检验')xlabel('时间')ylabel('呼出气计数率')%最小二乘法拟合曲线a=xlsread('C:\Users\谷柏辰\Desktop\data.xls');xdata=1:0.25:10;ydata=a(1,:);x0=[-5000 400 0.5]';figure(1);plot(xdata,ydata,'r*');hold on;[x,resnorm,residual]=lsqcurvefit('fun',x0,xdata,ydata); disp('系数矩阵x')xdisp('系数lamda=')x(3)disp('系数kesai=')kesai=x(1)*(x(3)-1.4808)/exp(9.1648)t=1:0.1:10;y=x(1)*exp(-1.4808*t)+x(2)*exp(-x(3)*t);plot(t,y)hold on;title('最小二乘法拟合');xlabel('时间');ylabel('头部计数率');t1=1;y1=x(1)*exp(-1.4808*t1)+x(2)*exp(-x(3)*t1)ynihe=x(1)*exp(-1.4808*xdata)+x(2)*exp(-x(3)*xdata) ychazhi=ynihe-ydata;figure(2);plot(xdata,ychazhi,'*');hold on;y0=0*xdata;plot(xdata,y0,'r','linewidth',2);title('差值分析');xlabel('时间');ylabel('拟合值与原始值求差');。