局部腦血流的測定數(shù)學建模

1、局部腦血流的測定一. 問題簡介 腦血流量是診斷和治療腦梗塞,腦出血,動脈瘤和先天性動,靜脈血管畸形等腦血管疾病的主要依據(jù)。測定腦血流量可為研究人腦在不同的病理和生理條件下的功能提供客觀指標，它對研究腦循環(huán)藥物的藥理作用也很有幫助。所以人們長期致力于尋找有效地測定腦血流量的方法。 近年來出現(xiàn)了以放射性同位素作示蹤劑測定人腦局部血流量的方法。這種方法大致可描述如下:由受試者吸入某種放射性同位素的氣體，然后將探測器置于受試者頭部某固定處，定時測量該處放射性同位素的計數(shù)率(簡稱計數(shù)率)，同時測量他呼出氣的計數(shù)率。 由于動脈血將肺部的放射性同位素輸送至大腦，使腦部同位素增加，而腦血流又將同位素帶離，使同

2、位素減少，實驗證明由腦血流引起局部地區(qū)計數(shù)率下降的速率與當時該處的記數(shù)率成正比，其比例系數(shù)反映了該處的腦血流量，被稱為腦血流量系數(shù)，只要確定該系數(shù)即可推算出腦血流量。動脈血從肺輸送同位素至大腦引起腦部計數(shù)率上升的速率與當時呼出氣的計數(shù)率成正比。 若某受試者的測試數(shù)據(jù)如下: 時 間(分) 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 頭部記數(shù)率 1534 1528 1468 1378 1272 1162 1052 947 848 757 674 599 呼出氣記數(shù)率 2231 1534 1054 724 498 342 2

3、35 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 頭部記數(shù)率 531 471 417 369 326 288 255 225 199 175 155 137 呼出氣記數(shù)率 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 頭部記數(shù)率 121 107 94 83 73 65 57 50 44 39 35 31 27 呼出氣記數(shù)率

4、 0 0 0 0 0 0 0 0 0 0 0 0 0 試建立確定腦血流系數(shù)的數(shù)學模型并計算上述受試者的腦血流系數(shù)。 備注:該題目是上海市(1990 年)大學生數(shù)學建模競賽A題。 二. 模型的假定1. 腦部計數(shù)率(記為ht()的上升只與肺部的放射性同位素有關(guān)，上升速度與呼出氣的記數(shù)率(記為 p()t )成正比，比例系數(shù)記為k ； 2. 腦部記數(shù)率ht()的下降只與該處腦血流量有關(guān)，其下降速度正比于ht()，比例系數(shù)為腦血流系數(shù)，記為K，這里忽略了放射性元素的衰變和其它因素； 3. 腦血流量在測定期間恒定，心臟博動，被測試者大腦活動，情感波動等帶來的變化忽略不予考慮； 4. 每次儀器測量為相互獨立

5、事件，各測量值無記憶相關(guān)； 5. 放射性同位素在人體內(nèi)傳遞是從吸入氣體(含有放射物)開始的，并假定一次吸入，因此認為同位素在肺中瞬時達到最大濃度； 6. 在吸入氣體瞬時，腦中放射物記數(shù)率為零； 7. 腦血流量與腦血流系數(shù)K成單值函數(shù)關(guān)系，求得后者即可確定前者。三. 模型的建立與分析 由于已知腦局部同位素的增減與已測定的頭部記數(shù)率ht()和呼出氣記數(shù)率 p()t 成正比關(guān)系,于是很自然地確定以腦部同位素量,即腦部記數(shù)率作為討論對象.。 1.原始模型的建立 設(shè)某時刻t 0時,頭部記數(shù)率為ht(),在t 時段后記數(shù)率ht()+t ，由假定可知, 頭部記數(shù)率的增量hht= ()(+t -ht)僅與三個

6、因素有關(guān)： (i) 肺動脈血將肺部的放射性同位素送至大腦,使腦部記數(shù)率增量為h ； (ii) 腦血流將同位素帶離,腦記數(shù)率下降為h ； (iii)放射性同位素自身有衰減引起記數(shù)率下降量為h ,設(shè)其半衰期為 . 因此,由醫(yī)學試驗及假定有 dh dh dh ln21 2 3 = kp(),t = Kh(),t =- ht(), dt dt dt 而 ht() = h ()t -h ()t +h ()t ， 123于是 dh dh dh dh123=-+ , dt dt dt dt即 dh ln2=-Kht()+kpt()- ht(). (dt 133 由于在測試時放射性同位素(如 Xe)的半衰期

7、一般很大,而測試時間又很短(大約十幾分鐘左右),由此假定 +,于是(1)式變?yōu)? dh=-Kh()t +kp()t . (dt 2. 算法模型的建立與改進 在建立算法模型之前,首先必須對 p()t 進行預測。作 p()tt 和-tln pt( ) t 的離散圖(圖1和圖2),由此發(fā)現(xiàn) p()t 與t 有近似于 Ae 的函數(shù)關(guān)系。 通過對ln pt( ) t 的離散圖2的觀察,去掉時刻6及6.5以后的樣本(這樣作的原因見文章后面的評注),再利用最小二乘法進行擬合得 lnpt( ) = 915158. -147577. t -147577. t其相關(guān)系數(shù)r = 0999887. ,由此得知 pt(

8、) = 9429.33e 3 -147577. t 我們作出 pt() = 9429.33e 的圖形,并將此圖和圖1放在一起圖3,由圖3及相關(guān)系數(shù)r = 0999887. 可以認為 p()t 確實是負指數(shù)曲線 -t pt() =Ae ,(A 9429.33, 147577. ). (2)及假設(shè)f，即h()0 = 0，解得 kA -tKtht() = (ee- ) ， (3) K - 此式從數(shù)學上來看并不復雜，但要利用此式求出參數(shù)K和k卻并非事，而參數(shù)K則需要在測試中使用，因此我們的問題歸結(jié)為：如何利實際測量值和(2)及h(0)=0去決定參數(shù)k和K這類數(shù)學問題稱為參數(shù)識問題。 下面建立幾個算法模

9、型: 算法模型.一般差分擬合法： 將方程(2)離散化，記時間步長為T，利用前插公式得： hh-nn+1 =-Kh +kp nnT hKThk= ()1- + Tp ， (4) nn+1 n中 hhtnTpptnT= (),()+ = + . nn002 用差分法求解，其截斷誤差為oT()，顯然大了些，為了提高精度準確度，最直接的方法是由插值方法獲得更多的結(jié)點，縮短步長，使斷誤差減少。如用三次樣條插值法在每兩個結(jié)點的中點進行插值，可2截斷誤差減少到原來的1/4，但仍然為oT()，且繼續(xù)縮短步長，計算量將成倍增加。 算法模型.改進的差分擬合法： 在這個算法中，我們注意方程(2)右端的線性項-Kh(

10、)t ，因此兩邊同Kt乘以e (積分因子)后可得: Ktde h()t Kt= kp()t e ， (5) dt對方程(5)利用差分離散化，并整理得： KTeh -hnn+1 = kp nT即： -KT KThehkTep=+ ， (6) n+1 n n-KT 2此時截斷誤差為oe( T )，顯然要比算法模型I 誤差要小，同時若將(6)-KT -KT中的e 展開，即eKT=-1 +oT()，略去高階無窮小，則得到： hKThk= ()1- + Tp nn+1 n這恰好是方程(4)，由此可見利用積分因子后得到了一個比模型I 精度要高的一個算法模型。 對于離散方程(4)或(6)可以通過聯(lián)立不同時刻

11、的方程組求得一系列K值，但是由于在實際測量中存在隨機誤差，以及離散化的截斷誤差,使得這些K值不盡相同。為了充分利用已測數(shù)據(jù)，我們利用最小二乘法擬合數(shù)據(jù)可得： hh= 0882488.+0078065p ， (7) nn+1 n在這里我們?nèi) =1，步長T = 025. ，擬合的復相關(guān)系數(shù)r = 09999997. 于0是將(7)與(4)式或(6)式比較可得參數(shù)K和k的值如下表所示： 算法模型 K的估計值 k的估計值 0.470048 0.31226 0.5004 0.353841 表1.算法算法的結(jié)果 上述兩個算法模型，計算簡單，但對誤差難以估計，并且對上述算法進行測試，兩個算法對K具有穩(wěn)定性

12、，而對另一個參數(shù)k卻不穩(wěn)定，同時也看到算法優(yōu)于算法，測試方法是預先假定一組K和k，按為未離散的公式(3)計算ht()在各時刻的值作為原始數(shù)據(jù)，再用差分公式和 最小二乘法求出K和k ，將它們與原假定值作比較，測試的結(jié)果見表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)式并作其連續(xù)圖,然后與離散圖作比較，同樣可以看出模型優(yōu)于模型(圖4對應于模型,圖5對應于模型)。 圖4 圖5 下面我們將給出另外一種算法對上述結(jié)果進

13、行改進. 算法摸型：線性迭代算法 如果設(shè)已給K和k的預測值K 和k ,記 0 0KK= + , 。 kk= + 0 0其中和稱為K和k相對于K 和k 的校正值(簡稱校正值)，將它們代0 0入(3)式并將右端關(guān)于和展開成Taylor 級數(shù)，同時略去和的二次及二次以上的項(即高階無窮小項)，得到 ()kA0 + -+tKt()0ht() = ()ee-K +-0 kA0 -tKt00A -tKt ()()ee-+ ee-+K - K -00-Kt t Kt00te ee-+Ak - 02K0 - ()K0 - = ht()利用理論值和實測數(shù)據(jù)誤差的平方和最小的原則來選取和,即選取和使 n 2 ()

14、hhtht=-() () iii=1最小.利用最小二乘法求得和后,較正K 和k 得 0 0 KK= + , kk= + 0 0將得到的新的參數(shù)K和k作為新的預測值，用同樣的方法繼續(xù)校正，直至和足夠小為止。 我們采用模型II的結(jié)果作為預測值,進行上述迭代程序得到的結(jié)論如表3所示: 迭代次 預測值K 預測值 校正值 校正值 ()h 數(shù) k 初始值 0.50004 0.353841 1 0.50004 0.353841 0.004573 0.066023 66.117 2 0.504613 0.419864 -0.000667 0.000025 65.3778 -6 -63 0.503946 0.4

15、19889 -1302. 10 -5459. 10 65.3557 -7 -74 0.503945 0.419884 -3053. 10 -4199. 10 65.3557 結(jié)果值 0.503945 0.419884 表3.算法的迭代結(jié)果 由表3可見算法模型的優(yōu)越性與準確性，并且得到K和k的最佳擬合值為： K=0.503945, k=0.419884 -7這種算法收斂速度很快，并且得到K值誤差數(shù)量級為10 四.模型的評價及注記 (1)我們所建立的前兩個算法模型計算簡單,但是穩(wěn)定性較差；第三個算法模型是穩(wěn)定的，并且具有快速收斂性，可獲得較精確的腦血流量系數(shù)K。利用得到的結(jié)果，(3)式和已知的數(shù)據(jù)作ht() t 的連續(xù)圖和離散圖，如圖