2020年第九屆數(shù)學(xué)建模國(guó)際賽小美賽

B題 血氧飽和度的變異性

原題再現(xiàn)：

??脈搏血氧飽和度是監(jiān)測(cè)患者血氧飽和度的常規(guī)方法。在連續(xù)監(jiān)測(cè)期間，我們希望能夠使用模型描述血氧飽和度的模式。
??我們有36名受試者的數(shù)據(jù)，每個(gè)受試者以1 Hz的頻率連續(xù)測(cè)試血氧飽和度約1小時(shí)。我們還記錄了參與者的以下信息，包括年齡、BMI、性別、吸煙史和/或當(dāng)前吸煙狀況，以及可能影響閱讀的任何重要疾病。
??我們想用這些數(shù)據(jù)來(lái)發(fā)現(xiàn)血氧飽和度變化的典型模式，這樣我們就可以用幾個(gè)參數(shù)來(lái)描述一個(gè)人。我們還想知道血氧飽和度序列的模式是否與年齡有關(guān)，即老年人與年輕人相比哪些特征發(fā)生了變化。理想情況下，這些特征應(yīng)具有生物學(xué)或醫(yī)學(xué)意義。

整體求解過(guò)程概述(摘要)

??脈搏血氧飽和度是監(jiān)測(cè)患者血氧飽和度的常規(guī)方法。脈搏血氧飽和度的使用有助于減少有創(chuàng)動(dòng)脈血?dú)夥治龊偷脱跹Y檢測(cè)的需要。一個(gè)可靠有效的血氧飽和度數(shù)學(xué)模型對(duì)進(jìn)一步研究具有重要意義。
??首先計(jì)算36名受試者血氧飽和度的均值和標(biāo)準(zhǔn)差，然后進(jìn)行初步的線性分析。結(jié)果表明，氧飽和度波動(dòng)較大，并伴有過(guò)飽和度。一般來(lái)說(shuō)，脈搏血氧飽和度顯示出很小的變異性。根據(jù)poincare曲線分析，血氧飽和度的變化主要由長(zhǎng)期變化組成。此外，我們還分析了平均SpO2與群體變異性之間的Pearson相關(guān)系數(shù)，發(fā)現(xiàn)SpO2水平與群體變異性呈負(fù)相關(guān)。DFA分析結(jié)果表明，時(shí)間序列是典型的分形時(shí)間序列，具有明顯的長(zhǎng)程相關(guān)性和長(zhǎng)程冪律。
??在此基礎(chǔ)上，探討了血樣多樣性各參數(shù)的具體模式，并提出利用ARMA時(shí)間序列模型對(duì)血氧飽和度的多樣性進(jìn)行建模和分析。我們對(duì)樣本進(jìn)行了單位根檢測(cè)，確定樣本為平穩(wěn)序列。然后對(duì)樣本的自相關(guān)函數(shù)（ACF）或偏自相關(guān)函數(shù)（PACF）進(jìn)行統(tǒng)計(jì)分析，確定模型的階數(shù)。最后通過(guò)機(jī)器學(xué)習(xí)得到ARMA模型的具體參數(shù)。通過(guò)殘差分析和D-W檢驗(yàn)驗(yàn)證了模型的正確性。通過(guò)模型分析可知，健康成人血氧飽和度濃度具有三階自相關(guān)和三階偏相關(guān)的特征。
??通過(guò)樣本熵分析、趨勢(shì)波動(dòng)分析和多尺度熵分析，探討了血氧飽和度序列模式與年齡的關(guān)系，得出以下結(jié)論：（1）年齡對(duì)平均血氧濃度無(wú)顯著影響。（2） 青年人和老年人的血氧飽和度變化是慢性的。（3） 從不同的尺度來(lái)看，老年人的樣本熵小于青年人，且在較高的尺度下差異更為明顯。從長(zhǎng)遠(yuǎn)來(lái)看，老化對(duì)OSV復(fù)雜性的降低有重要影響。
??綜上所述，該模型在血氧飽和度分析中準(zhǔn)確、真實(shí)，發(fā)現(xiàn)了年齡和血氧飽和度序列的具體特征，具有生物學(xué)或醫(yī)學(xué)意義。

模型假設(shè)：

??為了簡(jiǎn)化問(wèn)題并消除復(fù)雜性，我們做出以下假設(shè)。

??（1） 問(wèn)題中給出的數(shù)據(jù)是真實(shí)可靠的。該指令設(shè)置了一個(gè)限制，即提供的數(shù)據(jù)文件只包含我們應(yīng)該用于此問(wèn)題的數(shù)據(jù)，并且只有當(dāng)這些數(shù)據(jù)真實(shí)可靠時(shí)，我們的分析才有效。

??（2） 沒(méi)有其他影響因素。問(wèn)題中提供的數(shù)據(jù)涵蓋了可能影響研究人群OSV的所有重要醫(yī)療條件。

問(wèn)題重述：

??問(wèn)題背景
??脈搏血氧飽和度（pulseoximetry）是一種無(wú)創(chuàng)性測(cè)量血氧飽和度（SpO2）的技術(shù)。無(wú)論是在重癥監(jiān)護(hù)室、外科手術(shù)室，還是在一些門(mén)診，它都被證明是一種廣泛應(yīng)用的臨床方法。在這些環(huán)境中使用脈搏血氧飽和度有助于減少有創(chuàng)動(dòng)脈血?dú)夥治龊蜋z測(cè)低氧血癥的需要。
??利用血氧飽和度的變異性分析來(lái)進(jìn)一步測(cè)量血氧合的調(diào)節(jié)已引起越來(lái)越多的認(rèn)識(shí)。生理變異性分析的好處在于它可以為我們提供有關(guān)生理控制系統(tǒng)完整性的有用信息。氧飽和度變異性（OSV）分析可用于控制組織氧合監(jiān)測(cè)的心肺系統(tǒng)的完整性[13]。此外，它還用于睡眠呼吸紊亂的診斷，其中SpO2特征充分描述了SpO2調(diào)節(jié)，以識(shí)別睡眠呼吸紊亂的風(fēng)險(xiǎn)。在早產(chǎn)兒中，血氧飽和度變異性表現(xiàn)出明顯的特征，OSV穩(wěn)定增加，而平均SpO2值變化不大[8]。此外，研究人員試圖尋找OSV的診斷價(jià)值。例如，最近在孟加拉國(guó)一家三級(jí)醫(yī)院開(kāi)展的一項(xiàng)研究調(diào)查了實(shí)施OSV作為預(yù)測(cè)工具是否可以提高危重癥兒童的入院率。
??因此，一個(gè)穩(wěn)定有效的血氧飽和度數(shù)學(xué)模型將為進(jìn)一步的研究做出重要貢獻(xiàn)。

??問(wèn)題重述
???建立氧氣典型變化模式的數(shù)學(xué)模型，以確定人類(lèi)健康特征與OSV之間的關(guān)系。
???了解血氧飽和度序列的模式是否與年齡相關(guān)，以及與年輕人相比，老年人的哪些特征表現(xiàn)出明顯的變化。

模型的建立與求解整體論文縮略圖

全部論文請(qǐng)見(jiàn)下方“ 只會(huì)建模 QQ名片” 點(diǎn)擊QQ名片即可

部分程序代碼：(代碼和文檔not free)

li_ave=np.mean(li_np,axis=0) #Average of 36 people in 0-3500 seconds
li_ave_per=np.mean(li_np,axis=1) #Average blood oxygen per person
li_std_per=np.std(li_np,axis=1) #Blood oxygen variance per person
li_simp=li_np[0,:] #Data sample of the first person
#np.set_printoptions(threshold=np.inf)
print("li_std_per = ",li_std_per)
#Drawing
fig1 = plt.figure()
#Draw blood oxygen time series
plt.plot(range(len(li_simp)),li_simp)
plt.xlabel("Data points")
plt.ylim(90,105)
plt.ylabel("Oxygen saturation(\%)")
plt.title(’Oxygen Saturation Variability Over 1 Hour’)
#The relationship between mean blood oxygen and standard deviation
fig2 = plt.figure()
li_ave_std=np.std(li_ave_per)
li_ave_ave=np.mean(li_ave_per)
print("li_ave_ave = ",li_ave_ave)
print("li_ave_std = ",li_ave_std)
li_std_ave=np.mean(li_std_per)
li_std_std=np.std(li_std_per)
print("li_std_ave = ",li_std_ave)
print("li_std_std = ",li_std_std)
plt.scatter(li_ave_per,li_std_per)
plt.xlim(90,105)
plt.ylim(0,1.5)
plt.xlabel(’Mean SpO2 (\%)’)
plt.ylabel(’Standard Deviation of SpO2 (\%)’)
plt.title("Relationship between Mean Oxygen Saturation Level and Total Variability")
#Linear regression
model = LinearRegression()
model = model.fit(li_ave_per.reshape(-1,1), li_std_per)
plt.plot([93,100],[i*model.coef_+model.intercept_ for i in [93,100]])
#Correlation coefficient between mean blood oxygen and standard deviation
li_ave_std_r=np.mean(np.multiply((li_ave_per-np.mean(li_ave_per)),(li_std_per-np.mean(li_std_per))))/(np.std(li_std_per)*np.std(li_ave_per))
plt.text(92, 0.6, "r=\%.3f" \% li_ave_std_r)
fig3 = plt.figure()
for i in range(len(li_simp)-1):
plt.plot(li_simp[i], li_simp[i+1], color=’b’, marker=’o’)
plt.xlabel(’SpO2(n)(\%)’)
plt.ylabel(’SpO2(n+1)(\%)’)
plt.xlim(90,105)
plt.ylim(90,105)
plt.title("Poincare Plot for SpO2 data")
#SD calculation
all_SD1=[]
all_SD2=[]
for j in range(36):
SD1 = []
SD2 = []
li_temp=li_np[j,:]
for i in range(len(li_temp)-1):
SD1.append(li_temp[i+1]-li_temp[i])
SD2.append(li_temp[i+1]+li_temp[i])
ST1 = np.std(SD1)/np.sqrt(2)
ST2 = np.std(SD2)/np.sqrt(2)
all_SD1.append(ST1)
all_SD2.append(ST2)
SD1_ave=np.mean(all_SD1)
SD2_ave=np.mean(all_SD2)
SD1_std=np.std(all_SD1)
SD2_std=np.std(all_SD2)
print("SD1_ave = \%.2f"\%SD1_ave)
print("SD2_ave = \%.2f"\%SD2_ave)
print("SD1_std = \%.2f"\%SD1_std)
print("SD2_std = \%.2f"\%SD2_std)
plt.text(100,94,"SD1:\%.2f " \% all_SD1[0] + "\%")
plt.text(100,93,"SD2:\%.2f " \% all_SD2[0] + "\%")

function SampEnVal = SampEn(data, m, r)
data = data(:)’;
N = length(data);
Nkx1 = 0;
Nkx2 = 0;
for k = N - m:-1:1
x1(k, :) = data(k:k + m - 1);
x2(k, :) = data(k:k + m);
end
for k = N - m:-1:1
x1temprow = x1(k, :);
x1temp = ones(N - m, 1)*x1temprow;
dx1(k, :) = max(abs(x1temp - x1), [], 2)’;
Nkx1 = Nkx1 + (sum(dx1(k, :) < r) - 1)/(N - m - 1);
x2temprow = x2(k, :);
x2temp = ones(N - m, 1)*x2temprow;
dx2(k, :) = max(abs(x2temp - x2), [], 2)’;
Nkx2 = Nkx2 + (sum(dx2(k, :) < r) - 1)/(N - m - 1);
end
Bmx1 = Nkx1/(N - m);
Bmx2 = Nkx2/(N - m);
SampEnVal = -log(Bmx2/Bmx1);

function [mse,sf] = MSE_Costa2005(x,nSf,m,r)
% pre-allocate mse vector
mse = zeros([1 nSf]);
% coarse-grain and calculate sample entropy for each scale factor
for ii = 1 : nSf
% get filter weights
f = ones([1 ii]);
f = f/sum(f);
% get coarse-grained time series (i.e., average data within non-overlapping time
windows)
y = filter(f,1,x);
y = y(length(f):end);
y = y(1:length(f):end);
% calculate sample entropy
mse(ii) = SampleEntropy(y,m,r,0);
end
% get sacle factors
sf = 1 : nSf;

function F_n=DFA(DATA,win_length,order)
N=length(DATA);
n=floor(N/win_length);
N1=n*win_length;
y=zeros(N1,1);
Yn=zeros(N1,1);
fitcoef=zeros(n,order+1);
mean1=mean(DATA(1:N1));
for i=1:N1
y(i)=sum(DATA(1:i)-mean1);
end
y=y’;
for j=1:n
fitcoef(j,:)=polyfit(1:win_length,y(((j-1)*win_length+1):j*win_length),order);
end
for j=1:n
Yn(((j-1)*win_length+1):j*win_length)=polyval(fitcoef(j,:),1:win_length);
end
sum1=sum((y’-Yn).^2)/N1;
sum1=sqrt(sum1);
F_n=sum1;

全部論文請(qǐng)見(jiàn)下方“ 只會(huì)建模 QQ名片” 點(diǎn)擊QQ名片即可