2020年第九屆數學建模國際賽小美賽

B題 血氧飽和度的變異性

原題再現：

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

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

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

模型假設：

??為了簡化問題并消除復雜性，我們做出以下假設。

??（1） 問題中給出的數據是真實可靠的。該指令設置了一個限制，即提供的數據文件只包含我們應該用于此問題的數據，并且只有當這些數據真實可靠時，我們的分析才有效。

??（2） 沒有其他影響因素。問題中提供的數據涵蓋了可能影響研究人群OSV的所有重要醫療條件。

問題重述：

??問題背景
??脈搏血氧飽和度（pulseoximetry）是一種無創性測量血氧飽和度（SpO2）的技術。無論是在重癥監護室、外科手術室，還是在一些門診，它都被證明是一種廣泛應用的臨床方法。在這些環境中使用脈搏血氧飽和度有助于減少有創動脈血氣分析和檢測低氧血癥的需要。
??利用血氧飽和度的變異性分析來進一步測量血氧合的調節已引起越來越多的認識。生理變異性分析的好處在于它可以為我們提供有關生理控制系統完整性的有用信息。氧飽和度變異性（OSV）分析可用于控制組織氧合監測的心肺系統的完整性[13]。此外，它還用于睡眠呼吸紊亂的診斷，其中SpO2特征充分描述了SpO2調節，以識別睡眠呼吸紊亂的風險。在早產兒中，血氧飽和度變異性表現出明顯的特征，OSV穩定增加，而平均SpO2值變化不大[8]。此外，研究人員試圖尋找OSV的診斷價值。例如，最近在孟加拉國一家三級醫院開展的一項研究調查了實施OSV作為預測工具是否可以提高危重癥兒童的入院率。
??因此，一個穩定有效的血氧飽和度數學模型將為進一步的研究做出重要貢獻。

??問題重述
???建立氧氣典型變化模式的數學模型，以確定人類健康特征與OSV之間的關系。
???了解血氧飽和度序列的模式是否與年齡相關，以及與年輕人相比，老年人的哪些特征表現出明顯的變化。

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

全部論文請見下方“ 只會建模 QQ名片” 點擊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;

全部論文請見下方“ 只會建模 QQ名片” 點擊QQ名片即可