2022年亚太杯APMCM数学建模大赛A题结晶器熔剂熔融结晶过程序列图像特征提取及建模分析求解全过程文档及程序

2022年亚太杯APMCM数学建模大赛

A题 结晶器熔剂熔融结晶过程序列图像特征提取及建模分析

原题再现：

  连铸过程中的保护渣使钢水弯液面隔热，防止钢水在连铸过程中再次氧化，控制传热，为铸坯提供润滑，并吸收非金属夹杂物。结晶器熔剂的冶金性能主要取决于其在温度控制曲线下的熔化速率和结晶速率。因此，研究结晶器保护渣在结晶器壁和铸坯壳间隙中的相分布具有重要意义。

  连铸结晶器熔剂被添加到结晶器中液态钢的顶部。这些固体炉渣以粉末层的形式堆积在钢液表面，可以防止钢液温度下降过大而导致钢液液面结皮。然后，模具焊剂的温度逐渐上升到熔点，模具焊剂被熔化以形成烧结层。助熔剂的原料通过化学反应形成低熔点物质，然后形成液态渣，助熔剂的成分会发生一定的变化。这是一个熔化的过程。

  当保护渣完全熔化时，会在液态钢表面形成一层液态渣层。当液态炉渣从钢液面处的渣池渗入壳体与铜模壁之间的间隙时，就会形成渣膜。由于铸坯表面的高温，铸坯上的熔渣仍保持液相。然而，随着液态炉渣的温度在结晶器的纵向上随着铸坯表面的温度降低，渣膜在铜结晶器壁上被淬火固化，形成玻璃状固体渣膜（渣膜的固化行为），伴随结晶器的强制冷却，而渣膜在适当的条件下会在一定的区域结晶并形成一层晶层（渣膜的结晶行为），最终形成典型的三层渣膜结构：玻璃层、结晶层和液态渣层。这个过程就是结晶。

  由于高温、瞬态流体流动、复杂的相变和化学反应以及结晶器壁的不透明性，很难直接观察到结晶器熔剂的相变。SHTT-II型熔融结晶温度测试仪目前已广泛应用于模具熔剂结晶行为的观察。实验结束后，实验人员对图像进行逐一演示，将信息记录在图像左上角，并用肉眼和经验识别关键节点图像（见图1），以指导模具焊剂的设计，满足钢种的凝固要求。这一过程浪费了人力，阻碍了实验过程信息的开发。迫切需要开发序列图像的自动特征提取和数学建模技术。

  附件1中有562个模具助熔剂熔化和结晶的序列图像。这些序列图像是从实验开始时的第110秒到第671秒收集的。文件序列号遵循采集时间序列，每1s采集一次图像。附件1中的数字图像显示了这些信息（见图1）。每个图像的左上角标记有图像的相应时间以及1号热电偶和2号热电偶的温度值。

  为了实现模具熔剂熔化和结晶序列图像的自动特征提取和数学建模，请回答以下三个问题。

  问题1：使用图像分割和识别或其他技术，请自动提取每张图像左上角的1号热电偶和2号热电偶的温度，并将其自动导入附件2中相应的表格中（请编写一份循序渐进的技术操作文件），请制作温度-时间曲线图（1#线温-2#线温-时间图；1#线平均温度-2#线平均温度-时间图）。此外，1#线或2#线的测试结果不准确。请指出并解释。

  问题2：根据图1中的六个节点图像，应用数字图像处理技术，研究并量化模具熔剂熔化和结晶过程中相邻序列图像之间的动态差异。在此基础上，对量化后的不同特性进行时间序列建模，并根据数学模型的模拟结果讨论了结晶器熔剂的熔融结晶过程曲线。

  问题3：考虑到温度和时间的变化，以及问题2的研究结果，请建立一个数学模型来讨论温度和时间变化之间的函数关系以及熔化和结晶过程。

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

  通过建立像素判别模型，解决了图像中温度数据的识别问题。通过建立基于特征的状态方程模型，解决了图像特征的量化和模具熔融结晶过程曲线的建立问题。通过建立多项式回归模型，解决了温度和时间变化问题以及结晶器熔剂熔融与结晶过程的函数关系，并结合结晶动力学中的JMA公式验证了结果的可靠性。

  对于问题一，我们建立了像素判别模型，提取每张图片左上角的温度，并分别绘制1#和2#的温度-时间曲线。首先，对图像进行灰度处理，选择全局阈值对图像进行二值化，图像由0-1像素矩阵表示；然后，根据图片的数字字体和位置特征，将每个预处理后的图片分割成8个部分，然后通过滤波和降噪处理来提高识别率。然后，选择0-9数字和°C作为模板，通过计算图像和模板之间的相关系数，选择值最大的模板作为识别结果，并自动存储在表中；最后，使用卷积神经网络对结果进行了测试，结果是可靠的。

  对于问题二，我们建立了一个基于特征的状态方程模型，研究了熔融和结晶过程中的动态差异，并讨论了不同特征下熔融和结晶的过程曲线。首先，对图像进行预处理，将像素大小缩小到300*400；然后，使用“Canny”算子检测边缘，获得边缘轮廓图像，并使用find函数获得轮廓边缘图像的索引向量，以将模具与背景分离。然后，使用Tamura纹理特征算法提取粗糙度、对比度、方向性和线性，观察模具在熔化和结晶过程中的变化规律，并找到其关于时间的递归公式。最后，对平均灰度值和粗糙度进行可视化，实现多项式拟合，得到拟合后的熔融结晶过程曲线。通过粗糙度、对比度、方向性、线性和平均灰度值五个特征对结果进行测试和分析，证明了结果与实验现象的一致性，验证了基于特征的状态方程模型的可靠性。

  对于问题三，本文开发了一个多项式回归模型来解决将温度建立为时间的函数以及模具熔体的熔化和结晶的问题。首先，将多项式回归拟合到问题1的不同时刻的温度，以建立分段函数表达式。然后，基于问题2中发现的模具熔化和结晶的过程曲线，创建了在相同时刻熔化量与温度和结晶量与温度的多项式回归表达式。结晶速率与时间有关，熔融速率与时间相关。最后，结合结晶动力学，发现本文中发现的结晶速率与温度的关系图与JMA方程的图像高度一致，证明了所获得结果的准确性和可靠性。

模型假设：

  （1） 假设图像左上角的文本是顶部框架。

  （2） 假设图像记录的时刻是准确的，并且不考虑误差。

  （3） 假设模具焊剂正常，则使用过程没有问题。

问题分析：

  针对连铸过程中结晶器熔剂在高温加热下的复杂相变和化学反应，使直接观察结晶器熔剂的相态变得困难，提出了三个问题。首先，为了减少人力，有必要开发序列图像的数字识别和记录温度数据；然后，通过特征提取序列图像，并通过特征建立模型，研究模具熔融结晶过程曲线。最后，推导出提取的温度数据和模具熔融结晶过程曲线，推导出温度、熔融速率和结晶速率的函数关系表达式，确定了温度、熔融速度和结晶速率之间的函数关系式。这三个问题由浅入深，相互交织。

  问题一分析

  在对图像的分析中，在562个图像中的每个图像中要识别的数字位于图像的左上角，并且数据字体和数字之间的间隔是固定的。剪切将导致更小的图像大小和更少的像素点。设想图像将被放大10倍，为后续操作做准备。在剪切之前，本文设想对图像进行灰度处理和二值化，将图像简化为黑白双色，即两个值为0,1，将图像中的数字拆分，然后将数字剪切为大小为15893的像素，并通过均匀缩放将数字的大小设置为2828像素；本文注意到数字的边缘可能不平滑，并且为了使图像中数字的边缘更平滑；本文设想每个图像剪切出8个数字，但经观察，剪切的图像并非全部是数字，本文设想将未编号的°C识别为数字10，然后进行均匀处理。计算剪切图像数据和设计的标准数据之间的相关系数。根据预测，相关系数越大，图像数据与标准数据之间的相似性就越高。基于这一原理，本文设想建立一个像素判别模型，该模型将用于执行循环遍历，以识别每个图像与标准样本相比的温度，并导出，从而生成温度-时间分布图。本文设想对该模型进行比较分析，并使用卷积神经网络构建数字图像识别模型。如果得出的结论与像素判别模型的结论没有显著差异，那么像素判别模型是合理的，测试结果是准确的。

  问题二分析

  问题2是一种特征提取类型的问题。本文利用Tamura纹理特征提取了四个特征：粗糙度、对比度、方向性和线性。为了使计算更灵敏，图像采用灰度处理进行预处理，将图像变成黑白，然后将图像缩放为均匀大小。本文设想使用“Canny”算子来检测边缘，但边缘算子的局限性在于它只响应边缘并找到所有边缘，而非边缘应该被丢弃。因此，为了使边缘算子的结果更加准确，本文设想了一系列图像处理：首先对图像进行高斯平滑滤波和卷积，使图像边缘平滑；然后对梯度的幅度进行非最大值抑制，以使模糊的边界清晰；然后使用双阈值技术来进一步消除噪声。在处理后获得的图像中跟踪边缘，以获得边缘轮廓图像并将模具与背景分离。本文设想使用Tamura纹理特征提取，然后观察结晶和熔融过程中四个特征的变化，找出它们的变化模式，列出特征方程，并将计算后的结果绘制成图形，使结果更直观。

  问题三分析

  问题3属于函数关系派生类问题。基于问题1中的温度变化数据和问题2中的模具熔化和结晶的过程曲线，导出了它们之间的函数关系。为了函数表达式的准确性，我们选择使用多项式回归拟合，使离散数据与函数图像更相关。通过建立温度与时间、模具熔剂熔化与结晶与时间的函数关系，可以推导出模具熔剂熔化和结晶与温度的函数关系。将函数关系可视化并与结晶动力学公式相结合，以测试函数关系建立的准确性。如果图像相似，则结果是可靠的。

模型的建立与求解整体论文缩略图 (多种解法和文档)

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

程序代码：

部分程序如下:

clc
clear all
I=imread('D:\MATLAB\2022 亚太\Attachment 1\0110.bmp');
I1=rgb2gray(I);%灰度图像处理
I2=im2bw(I1,graythresh(I1));%二值图像处理
I3=I2(1:65,1:110);
I3=imresize(I3,10);
TT=I3(208:366,780:873);
ED_type='disk';
se=strel(ED_type,8);
TT=imopen(TT,se);
TT=imresize(TT,[28,28]);
TT=imcomplement(TT);
%%imwrite(TT,'a.jpg');
fulldate=dir('D:\MATLAB\2022 亚太\Attachment 1\*.bmp');
the_length=length(fulldate);
image0=fullfile('D:\MATLAB\2022 亚太\数字\new1#20572.bmp');
image0=imread(image0);
image1=fullfile('D:\MATLAB\2022 亚太\数字\new1#10132.bmp');
image1=imread(image1);
image2=fullfile('D:\MATLAB\2022 亚太\数字\new1#20177.bmp');
image2=imread(image2);
image3=fullfile('D:\MATLAB\2022 亚太\数字\new1#20189.bmp');
image3=imread(image3);
image4=fullfile('D:\MATLAB\2022 亚太\数字\new1#20203.bmp');
image4=imread(image4);
image5=fullfile('D:\MATLAB\2022 亚太\数字\new1#20251.bmp');
image5=imread(image5);
image6=fullfile('D:\MATLAB\2022 亚太\数字\new1#20594.bmp');
image6=imread(image6);
image7=fullfile('D:\MATLAB\2022 亚太\数字\new1#20589.bmp');
image7=imread(image7);
image8=fullfile('D:\MATLAB\2022 亚太\数字\new1#10645.bmp');
image8=imread(image8);
image9=fullfile('D:\MATLAB\2022 亚太\数字\new1#10110.bmp');
image9=imread(image9);
K{1}=image0;
K{2}=image1;
K{3}=image2;
K{4}=image3;
K{5}=image4;
K{6}=image5;
K{7}=image6;
K{8}=image7;
K{9}=image8;
K{10}=image9;
K{11}=TT;
%%比较相关系数
Z=[];
R=[];
for i = 1 : the_length
for j = 1:4
n=i+109;
image = strcat('D:\MATLAB\2022 亚太\数字
\new1#',num2str(j),'0',num2str(n),'.bmp');
Image = strcat('D:\MATLAB\2022 亚太\数字
\new2#',num2str(j),'0',num2str(n),'.bmp');
image = imread(image);
Image = imread(Image);
for k = 1 : 11
Z(k)=corr2(K{k},image);
Z2(k)=corr2(K{k},Image);
end
if Z(k)<=1
R(i,j)=find(Z==max(Z))-1;
else
R(i,j)=-1;
end
if Z2(k)<=1
R2(i,j)=find(Z2==max(Z2))-1;
else
R2(i,j)=-1;
end
end
end
%%找出没有图像的图片
wrong=[];
for i = 1 : the_length
if R(i,1)==-1
wrong=[wrong,i+109];
end
end
wrong = wrong';
for i = 1 : length(wrong)
for j = 1 : 4
image_path = strcat('D:\MATLAB\2022 亚太\数字
\','new1#',num2str(j),'0',num2str(wrong(i)),'.bmp');
image = imread(image_path);
for k = 1 : 11
Z(k)=corr2(K{k},image);
end
R(wrong(i)-109,j)=find(Z==max(Z))-1
end
end
%%热电偶温度变化
for i = 1 : the_length
if R(i,4)~=10
wendu_1(i)=R(i,1)*1000+R(i,2)*100+R(i,3)*10+R(i,4);
else
wendu_1(i)=R(i,1)*100+R(i,2)*10+R(i,3);
end
if R2(i,4)~=10
wendu_2(i)=R2(i,1)*1000+R2(i,2)*100+R2(i,3)*10+R2(i,4);
else
wendu_2(i)=R2(i,1)*1000+R2(i,2)*100+R2(i,3)*10+R2(i,4)-9;
end
end
wendu_1=wendu_1';
wendu_2=wendu_2';
%%校准
for i = 1 : the_length
if wendu_1(i)<300
wendu_1(i)=wendu_1(i)*10+1;
end
if wendu_2(i)<300
wendu_2(i)=wendu_2(i)*10+1
end
end
%%数据存入表 2
t=110:the_length+109;
t=t';
% xlswrite('Attachment 2.xlsx',t,'sheet1','B2:B563');
% xlswrite('Attachment 2.xlsx',wendu_1,'sheet1','C2:C563');
% xlswrite('Attachment 2.xlsx',wendu_2,'sheet1','D2:D563');
wendu_1=xlsread('Attachment 2.xlsx','C2:C563');
wendu_2=xlsread('Attachment 2.xlsx','D2:D563');
wendu_1(321,1)=1291;
%%画图
%%picture1
figure
%plot(t,wendu_1,'b',t,wendu_2,'r');
plot(t,wendu_1,'b','linewidth',1.5);
hold on;
plot(t,wendu_2,'r','linewidth',1.5);
grid on;
xlabel('Time/s');
ylabel('Temperature/℃');
title('1#wire temperature-2#wire temperature-time diagram');
%%picture2
pj_wendu_1=ones(length(wendu_1),1).*sum(wendu_1)/length(wendu_1);
pj_wendu_2=ones(length(wendu_2),1).*sum(wendu_2)/length(wendu_2);
figure
%plot(t,pj_wendu_1,'b',t,pj_wendu_2,'r');
plot(t,pj_wendu_1,'b','linewidth',1.5);
hold on;
plot(t,pj_wendu_2,'r','linewidth',1.5);
grid on;
xlabel('Time/s');
ylabel('Temperature/℃');
title('1#wire average temperature-2#wire average temperature-time
diagram');
%%compare
cnn_1=xlsread('D:\MATLAB\2022 亚太\温度.xls','B1:B562');
cnn_2=xlsread('D:\MATLAB\2022 亚太\温度.xls','C1:C562');
%%1#
figure
%plot(t,wendu_1,'b',t,cnn_1,'r');
plot(t,cnn_1,'k','linewidth',1.5);
hold on;
plot(t,wendu_1,'b','linewidth',1.5);
grid on;
xlabel('Time/s');
ylabel('Temperature/℃');
title('1#Forecast temperature comparison chart');
%%2#
figure
%plot(t,wendu_2,'b',t,cnn_2,'r');
plot(t,cnn_2,'k','linewidth',1.5);
hold on;
plot(t,wendu_2,'b','linewidth',1.5);
grid on;
xlabel('Time/s');
ylabel('Temperature/℃');
title('2#Forecast temperature comparison chart')；
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174

clc
clear all
close all
%%import data
digitDatasetPath = fullfile('E:\Matlab\matlab
apply\bin\shouxieshibie\shouxieshibie','/handwrite/');
imds = imageDatastore(digitDatasetPath,...
'IncludeSubfolders',true,'LabelSource','foldernames');%The folder
name is used as the data tag
%%Number of data set graphs
countEachLabel(imds)
numTrainFiles = 17;%Each number has 22 samples, and 17 samples are taken
as training data
[imdsTrain,imdsValidation] =
splitEachLabel(imds,numTrainFiles,'randomized');
%View the size of the image
img = readimage(imds,1);
size(img)
%%Define the structure of convolutional neural network
layers = [
%Layer of input
imageInputLayer([28 28 1])
%Layer of convolution
convolution2dLayer(5,6,'Padding',2)
batchNormalizationLayer
reluLayer
maxPooling2dLayer(2,'stride',2)
convolution2dLayer(5,16)
batchNormalizationLayer
reluLayer
maxPooling2dLayer(2,'stride',2)
convolution2dLayer(5,120)
batchNormalizationLayer
reluLayer
%Final layer
fullyConnectedLayer(11)
softmaxLayer
classificationLayer];
%%Trained neural network
%%Setting training Parameters
options = trainingOptions('sgdm',...
'MaxEpochs',50,...
'ValidationData',imdsValidation,...
'ValidationFrequency',5,...
'Verbose',false,...
'Plots','training-progress');%Shows the progress of the training
%Train the neural network, save the network
net = trainNetwork(imdsTrain,layers,options);
save 'CNNet.mat' net
%load('CNNet.mat')
%%Tagged data
mineSet = imageDatastore('E:\Matlab\matlab apply\bin\手写识别\手写识别
\hw22\hw22', 'IncludeSubfolders',false,'FileExtensions','.jpg' );
mLabels = cell(size(mineSet.Files,1),1);
for i = 1:size(mineSet.Files,1)
[filepath,name,ext] = fileparts(char(mineSet.Files{i}));
mLabels{i,1} = char(name);
end
mLabels2 = categorical(mLabels);
mineSet.Labels = mLabels2;
%%Use the network to classify and calculate accuracy
%Handwritten data
YPred = classify(net,mineSet);
YValidation = mineSet.Labels;
%Rate of accuracy of calculation
accuracy = sum(YPred == YValidation)/numel(YValidation)
%Plot the forecast results
figure;
nSample = 11;
ind = randperm(size(YPred,1),nSample);
for i = 1:nSample
subplot(2,fix((nSample+1)/2),i)
imshow(char(mineSet.Files(ind(i))))
title(['Forecast：' char(YPred(ind(i)))])
if char(YPred(ind(i))) ==char(YValidation(ind(i)))
xlabel(['True:' char(YValidation(ind(i)))],'Color','b')
else
xlabel(['True:' char(YValidation(ind(i)))],'color','r')
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

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