2022年华中杯数学建模挑战赛C题矿井提升机钢丝绳的缺陷分析求解全过程文档及程序

2022年华中杯数学建模

C 题 矿井提升机钢丝绳的缺陷分析

原题再现：

  矿井提升机是一种大型提升机械设备。由电机带动机械设备，以带动钢丝绳从而带动容器在井筒中升降，完成输送任务。按钢丝绳在卷筒上的连接形式分为缠绕式提升机、摩擦式提升机，按绳的多少可分为二绳、四绳、六绳提升机。图 1 为矿井提升机，图 2 为四绳提升机现场应用场景图。

  目前，基于“磁平衡”的检测技术被广泛地应用到提升机钢丝绳的全面安全监控中，该技术的优势就是在不损害钢丝绳性能及不影响钢丝绳使用的前提下，通过传感器获取监控数据信号，并通过数据信号进行分析，监测钢丝绳中是否出现断丝、变形、磨损、锈蚀等缺陷，同时给出缺陷的位置、大小（程度）、性质（类型）和发生缺陷的数量等信息，据此判定钢丝绳所处的技术状态（如合格与否、寿命等），进而达到对钢丝绳的全面自动在线安全监控，提高维护效率，节约维护成本。图 3 给出了理想状态下的监测数据的一段波形图，并标注了 5 种缺陷对应的波形，也给出了 5 种缺陷的量化定义，即为电压值和平均值之差的绝对值与平均值的百分比。

  附件 1 给出了某六绳提升机六根长度为 960 米的钢丝绳在井筒中升降往返 10 次的监测数据，数据相关说明见附件 2。请依据附件 1 数据，回答以下问题：
  问题 1：由于钢丝绳有缺陷（如断丝、变形、磨损、锈蚀等），会导致监测数据产生波动，又监测时钢丝绳会有某种程度的晃动等因素的影响，监测数据会出现不同程度的噪音，导致误报。对于附件 1 给出的监测数据，请建立数学模型及设计算法去除噪音，并找出六根钢丝绳经过 10 次往返运行后所有缺陷的数量、程度及位置。
  问题 2：结合问题 1 的结果，建立数学模型，给出钢丝绳安全性能的量化评价方法，依据附件 1 数据，计算六根钢丝绳经过 10 次往返运行之后的安全性能评价值。

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

  钢丝绳作为矿井提升系统牵引、运输、提升过程中的关键部件，在其长期工作过程中会受到矿井复杂环境的腐蚀、磨损、物理撞击、挤压等各种因素的影响，会出现断丝、磨损、变形、断股等不同程度的损伤和缺陷，导致使用强度和寿命降低，使得钢丝绳发生骤断的风险升高，威胁到人员财产生命安全，造成生产停顿。因此，对钢丝绳具体的损伤和缺陷进行定量识别，对于评价钢丝绳安全性能和保障生产安全具有重要的理论意义和实用价值。
  在问题一中，我们首先计算出上行和下行监测数据的每一帧的具体位置，接着对 10 次往返数据进行了数据整合，把在同一位置测得的不同数据进行平均处理，并进行重新排序。考虑到实际监测数据的复杂性，我们引入了多尺度小波分解滤波的思想，并利用MATLAB小波分析工具箱对钢丝绳检测信号进行分析。通过对原始信号进行 12 尺度分解，设置各分解信号的阈值完成抑噪处理；然后对信号进行重构，得到了去噪后的信号数据，在此基础之上，计算得到了 6 根钢丝绳的缺陷数目、程度和位置数据。
  在问题二中，我们根据已经得到的 6 根钢丝绳的缺陷数据，把影响钢丝绳安全性能指标的主要因素分为缺陷数目、缺陷程度（微小断丝、变形等）、缺陷位置（密集、稀疏等）和电压值数据的稳定性四类。从具体的数据出发，分别建立了相应的量化评价模型，给出了各钢丝绳各因素的评分指标。之后，对各钢丝绳各因素的评分指标应用熵权法模型，得到了最终的钢丝绳安全性能评分；结果显示，1 号钢丝绳的安全性能非常差，2 号钢丝绳的安全性能相对较差，3 号钢丝绳的安全性能相对较好，4、5 和 6 号钢丝绳的安全性能非常好。另外，我们又通过 TOPSIS 方法对各钢丝绳各因素的评分指标进行处理，其结果显示的相对接近度与安全性能评分十分接近，互相验证了两种模型的有效性。

问题分析：

  在对问题的分析中，考虑到实际工程的安全性，要精确的确定六根钢丝绳有缺陷的位置和程度，而附录中只给出了在对应脉冲下检测到钢丝绳电压，没有准确地定位到钢丝绳的具体位置；所以需要先对钢丝绳的数据进行处理，来确定钢丝绳每一帧数所对应的具体位置。另外，10 次的往返测量必然会导致在相同位置的重复测量，可能导致这些位置对应有不同的数据，为了减小误差和数据冗余，有必要对重复测量的数据进行处理。
  对于问题一，我们考虑到往返测量中监测数据不是从同一端记录的情况，分别计算出了上行和下行所测出的每一帧数据所对应的具体位置。接着把五次往返的数据（上行和下行的监测数据）进行整合，把在同一位置的不同数据进行平均处理，之后按照位置大小升序排列。由于噪音成分会对监测的数据造成干扰，我们采用对非平稳信号具有良好处理性能离散小波变换进行除噪。考虑到实际监测数据的复杂性，我们引入了多尺度小波分解滤波的思想，并通过对 MATLAB 小波分析工具箱的合理运用，得到了除噪后的监测数据。最后，我们根据 5 种缺陷量化的定义，设计算法找出了六根钢丝绳经过 10 次往返运行后所有缺陷的数量、程度和位置。
  对于问题二，在问题一中得到的结果的基础上，我们把影响钢丝绳安全性能的指标大致分为以下几类，微小断丝缺陷数目，变形缺陷数目，轻度断丝数目，中度断丝数目，内部断股数目，各类型缺陷所在的位置，任意钢丝绳所测得电压值的稳定性。通过定义指标，归一化处理，加权平均来得到相应的评分，之后通过熵权法计算综合的评价指数。

问题一模型的建立与求解

  在检测过程中，受到矿井复杂恶劣环境的影响，探头传感器输出的信号中掺杂了很多背景噪声信号，其中，断丝缺陷信号是一随机出现的局部异常信号，这是钢丝绳断丝缺陷定量识别的主要依据。所以，消除背景噪声信号是对钢丝绳断丝缺陷定量识别的前提。为确定断丝缺陷的位置，需要保持检测信号在空间域上的识别性，同时对其各个频段上进行抑噪处理。而传统的傅里叶变换在频域的分辨率为无穷大，在时域的分辨率为零，使得傅氏分析对损伤信号的处理、准确定位比较困难。为了克服傅氏分析的缺点，我们采用小波变换的方法对信号进行处理，小波变换是一个“时间—频率”窗口，窗口的形状随频率改变，但大小不变，是进行信号时频分析和处理时频局域化的分析方法。相比于傅里叶变换，小波变换保留了时间和频率的分辨率，更适合处理钢丝绳损伤这种突变的信号。离散小波变换是对基本小波的尺度和平移进行离散化，在处理非平稳信号时具有更好的性能。下面简要介绍离散小波变化的基本理论。

  小波滤波的思想是将外部断丝信号进行离散小波分解，得到低频系数和高频系数，噪声一般包含在高频系数中，所以对高频系数做阈值量化处理，再进行系数重构，得到去除噪声后的信号。实际中的检测信号较为复杂，单次小波分解不能很好的滤除噪声信号，通常采用多尺度小波分解滤波的方法，对检测信号进行多次分解。具体结构图如下所示：

  在前面我们简要探讨了多尺度小波分解滤波的方法。在滤波的过程中，噪声的高频部分随着尺度的增大迅速减小，但缺陷信号的奇异部分在小波变换大尺度上信号被保留[3]。小波分解的层次越多，意味着分解的尺度越大，对高频部分的滤除效果就会越好；但如果小波分解的尺度过大，就有可能造成部分有用信号的损失，且分解所用的时间过长，效率大大降低。因此，选择一个合适的小波分解尺度是进行有效滤波的关键。
  而在实际的工程中，监测信号是复杂多变的，难以快速地选择小波滤波分解的合适尺度；所以小波滤波分解尺度的选择是要结合实际的监测数据，经过多次实验来确定的，采用编程的方式来实现小波滤波的分解具有一定的复杂性，且实际操作不便。MATLAB 是一款强大的数据分析处理软件，为我们提供了小波分析的专用工具箱；本着“简单易行”的实际工程处理原则，为了得到更好的滤波效果，考虑采用 MATLAB 中的工具箱进行反复的实验，以得到最为优质的滤波信号。下面以 1 号钢丝绳为例，给出 MATLAB 工具箱所分析的结果。对 1 号钢丝绳的电压值信号进行 12 层小波分解，结果如图：

  经过 12 层小波分解，可以完全将损伤漏磁信号从各种信号中分离出来。高频干扰主要集中在信号 d1、d2、d3 、d4 上，在小波变换中对应着小的尺度，其能量随尺度的增大迅速减小，高频干扰主要是由钢丝绳股波和大功率电机引起；而低频扰动主要分布在 a12 上，其主要是由钢丝绳摆动等因素造成；损伤信号则主要集中 d9、d10、d11、d12 上。通过对 d1~d12 的数据进行阈值设置，得到去噪后的信号波形如图所示：

  通过原始信号与重构信号对比图可以明显地看出，小波分析方法在消噪的效果非常好，去除掉了大部分的干扰噪声，同时保留了大部分的缺陷信号。从图中幅值的变化和大小可以定性的看出，一号钢丝绳的缺陷最多，而且最严重，2 号、3 号和 4 号钢丝绳的缺陷相对较少，5 号、6 号钢丝绳几乎没有缺陷。利用 MATLAB 导出的重构信号的电压值数据，总的数据量为 2942635 帧，计算得到的各钢丝绳的缺陷程度及对应数目如下表所示：

问题二模型的建立与求解

  监测钢丝绳的目的就是为了更好地了解钢丝的实时性能，保障矿井提供机工作的安全。在前面我们对钢丝绳的监测数据进行了处理，并找到了六根钢丝绳经过 10 次往返运行后所有缺陷的数量、位置和程度。为了更好地了解钢丝绳的安全情况，我们把影响钢丝绳的安全性能指标主要分为以下几种，即微小断丝缺陷数目，变形缺陷数目，轻度断丝数目，中度断丝数目，内部断股数目，各类型缺陷所在的位置，任意钢丝绳所测得电压值的稳定性等。下面将基于获得的数据对这些安全性能指标进行量化分析，并对钢丝绳的安全性能做出评价。首先，我们考虑钢丝绳缺陷的数量和程度对钢丝绳安全性能的影响，我们知道钢丝绳缺陷的程度越严重，对实际工程安全性的威胁就越大，也就是缺陷程度越严重的钢丝绳的危险程度越高；在实际中，钢丝绳的危险程度和缺陷程度之间的关系往往是非线性的，所以我们考虑用指数函数去刻画这种非线性关系：

论文缩略图：

程序代码：

load SAA.mat
load LLL.mat
%将缺陷程度分别定义为1~5级
c1=[];%用于存放缺陷位置
c2=[];
c3=[];
c4=[];
c5=[];
d1=[];
d2=[];
d3=[];
d4=[];
d5=[];
k=[];
%设处理后有n个位置
%a(i)为处理后的电压
%b(i)为处理后的位置
a=SAA;
s=0;
n=length(a(:));
b=LLL;
for i=1:n
 %初始化一个数组,用于传递值
s=a(i)+s; %计算电压之和
end
a0=s/n; %计算电压均值
num1=0;
num2=0;
num3=0;
num4=0;
num5=0;%用于存放缺陷数量
for i=1:n 
 k(i)=(abs(a(i)-a0))/a0;
 if(k(i)<0.01)
 continue; %没有缺陷，继续循环
 elseif(k(i)>=0.01 && k(i)<0.016) 
 num1=num1+1;
 c1(end+1)=b(i);
 d1(end+1)=a(i);
 elseif(k(i)>=0.016 && k(i)<0.025)
 c2(end+1)=b(i);
 num2=num2+1;
 d2(end+1)=a(i);
 elseif(k(i)>=0.025 && k(i)<0.05)
 c3(end+1)=b(i);
 num3=num3+1;
 d3(end+1)=a(i);
 elseif(k(i)>=0.05 && k(i)<0.14)
 c4(end+1)=b(i);
 num4=num4+1;
 d4(end+1)=a(i);
 else
5(end+1)=b(i);
 num5=num5+1;
 d5(end+1)=a(i);
 end
 NUM=num1+num2+num3+num4+num5; 
end
plot(b,a,'.',c1,d1,'s',c2,d2,'+',c3,d3,'^',c4,d4,':',c5,d5,'*')

load CAA.mat
load CLLL.mat
a=CAA;
%将缺陷程度分别定义为 1~5 级
c1=[];%用于存放缺陷位置
c2=[];
c3=[];
c4=[];
c5=[];
d1=[];
d2=[];
d3=[];
d4=[];
d5=[];
k=[];
%设处理后有 n 个位置
%a(i)为处理后的电压
%b(i)为处理后的位置
s=0;
n=length(a(:));
b=CLLL;
for i=1:n
 %初始化一个数组,用于传递值
s=a(i)+s; %计算电压之和
end
a0=s/n; %计算电压均值
num1=0;
num2=0;
num3=0;
num4=0;
num5=0;%用于存放缺陷数量
for i=1:n 
 k(i)=(abs(a(i)-a0))/a0;
 if(k(i)<0.01)
 continue; %没有缺陷，继续循环
 elseif(k(i)>=0.01 && k(i)<0.016) 
 num1=num1+1;
 c1(end+1)=b(i);
 d1(end+1)=a(i);
 elseif(k(i)>=0.016 && k(i)<0.025)
 c2(end+1)=b(i);
 num2=num2+1;
 d2(end+1)=a(i);
 elseif(k(i)>=0.025 && k(i)<0.05)
 c3(end+1)=b(i);
 num3=num3+1;
 d3(end+1)=a(i);
 elseif(k(i)>=0.05 && k(i)<0.14)
 c4(end+1)=b(i);
 num4=num4+1;
 d4(end+1)=a(i);
 else
 c5(end+1)=b(i);
 num5=num5+1;
 d5(end+1)=a(i);
 end
 NUM=num1+num2+num3+num4+num5; 
end
plot(b,a,'.',c1,d1,'s',c2,d2,'+',c3,d3,'^',c4,d4,':',c5,d5,'*')
xlabel('位置/m')
ylabel('电压值/v')
legend('整个钢丝绳','微小断丝','变形','轻度断丝','中度断丝','内部断股')
l1=0;
for o1=1:num1
 l1=l1+(d1(o1)-480).^2;
 L1=l1;
end
l2=0;
for o2=1:num2
 l2=l2+(d2(o2)-480).^2;
 L2=l2;
end
l3=0;
for o3=1:num3
 l3=l3+(d3(o3)-480).^2;
 L3=l3;
end
l4=0;
for o4=1:num4
 l4=l4+(d4(o4)-480).^2;
 L4=l4;
end
l5=0;
for o5=5:num5
 l5=l5+(d5(o5)-480).^2;
 L5=l5;
end
LL=log(L1+L2+L3+L4+L5);
var(a,1)

获取论文及程序见下方