2023年认证杯SPSSPRO杯数学建模C题(第一阶段)心脏危险事件全过程文档及程序

2023年认证杯SPSSPRO杯数学建模

C题 心脏危险事件

原题再现：

  心脏的每一次搏动都伴随着心脏的电生理活动。心脏的起博点通过放电，使电流传导到每个心肌纤维，接收到电信号后，相应的心肌纤维完成一次收缩，心脏也就随之搏动一次。而心脏的电信号可以传导到体表皮肤，并且不同体表部位所检测到电信号表现不同。这样，在体表的特定部位放置电极，通过心电图机，可以记录到心电数据。对患有严重心脏疾病的人来说，心电的实时监测是检测心律失常的重要手段。
  为使心电监测更加有效，心电图机应当在心电图产生异常时能够做到实时报警。所以我们需要在很短时间内对心律失常进行正确的判断。我们在已有的心电图数据中找到了一些有代表性的片段，其中有正常心搏，也有多种心律失常的情况。每个片段长度为 2 秒。在数据文件中，我们记录的是心电波形的功率谱密度，从 0 Hz 到 180 Hz，频率间隔为 0.5 Hz。也就是第一行记录的是 0 Hz（直流分量）的数据，第二行记录的是 0.5 Hz，第三行记录的是 1Hz，依此类推。
  第一阶段问题：
  1. 请你和你的团队建立有效的数学模型，将所给的数据文件进行分类。除正常心搏外，请将心律失常的情况分为不同的类别，并指明类别的总数；
  2. 请给出每种心律失常类型的判断标准，以便我们能够核实判断方法的生理学意义，并将判断方法应用到临床监测设备上；
  3. 某些类型的心律失常一旦发生，心脏立即失去供血功能，此时病人的情况极为危急。另外一些类型的心律失常则不会如此危险，我们可以有稍多一些的救治时间。请参考正常的心搏过程，估计每种心律失常情况的危险程度，按照危险程度对数据文件进行粗略的排序或分级。

整体求解过程概述(摘要)

  本文首先通过对附件中存放的 869 组心律失常以及 147 组心律正常的心电波形的功率谱密度数据进行分析，利用 K-means 和 DBSCAN 两种聚类算法对心律失常数据进行了聚类划分，最终确定为 4 类心电功率谱。之后我们给出了心律失常数据聚类划分依据的特征参数的概念以及其对应的生理意义和临床诊断价值。最后，将心律失常数据与正常心律数据进行对比分析，将不同心律失常情况按照规定的标准进行了危险等级的简单划分。
  针对问题一，我们首先将数据导入到 Matlab 中，分析心律失常的心电波形的功率谱密度数据，对所有数据进行集中观察，发现 7 条数据的直流分量及主峰异常大，将其作为后续分类是否有效的参考。之后我们提取了基频、谱峰大小、谱峰数量、平均值、标准差等 9 个特征参数并进行了归一化及标准化的预处理，利用 K-means 和 DBSCAN 两种聚类算法对心律失常数据进行了聚类划分并进行结果对比及分析。
  针对问题二，我们根据所选特征参数的概念以及功率谱密度的意义和心脏活动的过程，对各特征参数进行了介绍以及其对于临床检测意义的说明。将功率谱密度与常规时域心电图判断心律失常类型的方法进行对比研究，并判断所划分的几种心律失常类型的可能病症。
  针对问题三，我们首先将正常数据和异常数据进行同样的特征参数提取以及预处理，将众多正常数据总结提取出一条最具代表性的数据，将异常数据与之进行差异性计算，根据差异结果进行心律异常数据的危险等级划分和排序。

问题分析：

  依照心脏起搏过程，心电图机可以很好地记录心电信号。对于患有严重心跳的病人，心电实时检测是很重要的。为在短时间做到心电信号的检测并报警，可以进一步观察短时间电信号片段的功率谱密度。在研究过程中，我们将充分利用现有心电图数据中的有代表性片段，包括正常心搏和多种心律失常情况。每个片段长度为 2 秒，数据文件记录的是心电波形的功率谱密度，从 0 Hz 到 180 Hz，频率间隔为 0.5 Hz。我们将依据这些数据，结合先进的计算机技术，开发出自动识别心律失常类型的模型，为心电监测提供实时报警功能，以提高诊断效率和准确率。
  问题一：将所给心律异常数据进行分类。
  思路分析：首先将数据导入到 Matlab 中，分析心律失常的心电波形的功率谱密度数据，将所有数据画在同一张图中进行集中观察。根据观察结果提取分类所需特征参数并进行一定的预处理，之后利用聚类算法对心律失常数据进行了聚类划分并进行结果对比分析。
  问题二：给出每种心律失常类型的判断标准，以便核实判断方法的生理学意义并实际应用。
  思路分析：根据所选特征参数的概念，结合功率谱密度的意义和心脏活动的过程对各特征参数进行介绍，说明对于临床检测的意义。将功率谱密度图与常规时域心电图判断心律失常类型的方法进行对比研究，判断所划分的几种心律失常类型的可能病症。若要应用到实际检测中，可将分类结果作为数据集利用 BP、CNN 等神经网络等进行模型训练，用于对实际病症的辅助诊断。
  问题三：参考正常的心搏过程，估计每种心律失常情况的危险程度，并按照危险程度对数据文件进行粗略的排序或分级。这将有助于医生迅速判断患者的病情严重程度，并采取相应的救治措施。
  思路分析：将心律正常数据进行和异常数据同样的特征参数提取以及预处理，将处理后的众多正常数据总结提取出一条最具代表性的数据，将异常数据与之一一进行差异计算，根据差异结果进行心律异常数据的危险等级划分和排序。

模型假设：

  1.假设题目所给的数据真实可靠
  2.假设数据具有普遍性和适用性
  3.假设心电图功率谱密度可以完全反映心率失常的情况
  4.假设病人心电图检查时未受外界的影响
  5.假设模型可以反映所有心律失常情况

完整论文缩略图

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：(代码和文档notfree)

clc
close all
clear
% X1=load('C:\All1\2023.csv');
% X2=load('C:\All1\20234.csv');
C1=importdata('C:\All1\20345.txt');
C2=importdata('C:\ist.txt');
% C1(1)
for c1=1:length(C1)
 X1(:,c1)=load(strcat('C:\thm\',cell2mat(C1(c1))));
end
for c2=1:length(C2)
 X2(:,c2)=load(strcat('C:\bnormal\',cell2mat(C2(c2))));
end
% for i1=1:length(C1)
% X1_max=max(X1(2:end,i1));
% X1(:,i1)=X1(:,i1)/X1_max;
% end
% for i2=1:length(C2)
% X2_max=max(X2(2:end,i2));
% X2(:,i2)=X2(:,i2)/X2_max;
% end
figure
for i=1:length(C1)
% subplot(211)
 plot(0:0.5:180,10*log(X1(1:end,i)),'LineWidth',1);
 xlabel('频率(Hz)')
 ylabel('物理功率谱密度(dBW)')
 set(gca,'FontSize',14,'FontName','Times New Roman')
 grid on
 subplot(212)
 plot(X2(2:25,i));
end
X11=X1(2:100,:);
X22=X2(2:100,:);
for ii=1:length(X11)
 plot(X1(2:100,ii));
 [pks1,loc1]=findpeaks(X11(:,ii));
 num_peak(ii)=length(loc1);
 for i=1:length(pks1)-1
 range(i)=loc1(i+1)-loc1(i);
 end
 range_peak(ii)=mean(range);
 loc_max(ii)=loc1(pks1==max(pks1));
 pks_max(ii)=max(pks1);
 base_fre(ii)=loc1(1)*0.5;
 energy(ii)=sum(X1(2:100,ii));
 mean_cs(ii)=mean(X1(2:100,ii));
 %stan_cs(ii)=sqrt(sum(X1(2:25,ii)-mean_cs(ii))/(25-2));
 DC(ii)=X1(1,ii);
 stan_cs(ii)=std(X1(2:100,ii));
 
Feature(:,ii)=[num_peak(ii);range_peak(ii);loc_max(ii);pks_max(ii);base_fre(ii);energy(ii);DC(ii);stan_cs(ii)
];
 
end
Feature=[Feature;ones(1,147)];
for ii2=1:length(X22)
plot(X2(2:100,ii2));
 [pks2,loc2]=findpeaks(X22(:,ii2));
 num_peak2(ii2)=length(loc2);
 for i=1:length(pks2)-1
 range2(i)=loc2(i+1)-loc2(i);
 end
 range_peak2(ii2)=mean(range2);
 loc_max2(ii2)=loc2(pks2==max(pks2));
 pks_max2(ii2)=max(pks2);
 base_fre2(ii2)=loc2(1)*0.5;
 energy2(ii2)=sum(X2(2:100,ii2));
 DC2(ii2)=X2(1,ii2);
 stan_cs2(ii2)=std(X2(2:100,ii2));
 ratio_max2(ii2)=X22(loc2(pks2==max(pks2))+1)/max(pks2);
 %Feature2(:,ii2)=[num_peak2(ii2);range_peak2(ii2);loc_max2(ii2);pks_max2(ii2);base_fre2(ii2);energy2(ii2);
stan_cs2(ii2)];
end
max1=max(num_peak2);
num_peak2=num_peak2/max1;
max2=max(range_peak2);
range_peak2=range_peak2/max2;
max3=max(loc_max2);
loc_max2=loc_max2/max3;
max4=max(pks_max2);
pks_max2=pks_max2/max4;
max5=max(base_fre2);
base_fre2=base_fre2/max5;
max6=max(energy2);
energy2=energy2/max6;
max7=max(stan_cs2);
stan_cs2=stan_cs2/max7;
max8=max(DC2);
DC2=DC2/max8;
max9=max(ratio_max2);
ratio_max2=ratio_max2/max9;
Feature2=[num_peak2;range_peak2;loc_max2;pks_max2;base_fre2;energy2;stan_cs2;DC2;ratio_max2];
%stan=[std(num_peak2),std(range_peak2),std(loc_max2),std(pks_max2),std(base_fre2),std(energy2),std(stan_cs2),
std(DC2),std(ratio_max2)];
%D= pdist(Feature2, 'seuclidean',stan);
%Feature2=[Feature2;zeros(1,869)];
%feature=[Feature,Feature2];
% figure
% for i=1:length(C1)
% X11=0.5*[flip(X1(2:end,i));X1(:,i)];
% plot(X11)
% plot(abs(ifft(X11)))
% 
% end

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可