1 题目简介

题目来源于2019 高教社杯全国大学生数学建模竞赛D题——空气质量数据的校准。
空气污染对生态环境和人类健康危害巨大，通过对“两尘四气”（PM2.5、PM10、CO、NO2、SO2、O3）浓度的实时监测可以及时掌握空气质量，对污染源采取相应措施。虽然国家监测控制站点（国控点）对“两尘四气”有监测数据，且较为准确，但因为国控点的布控较少，数据发布时间滞后较长且花费较大，无法给出实时空气质量的监测和预报。某公司自主研发的微型空气质量检测仪（如图所示）花费小，可对某一地区空气质量进行实时网格化监控，并同时监测温度、湿度、风速、气压、降水等气象参数。

由于所使用的电化学气体传感器在长时间使用后会产生一定的零点漂移和量程漂移，非常规气态污染物（气）浓度变化对传感器存在交叉干扰，以及天气因素对传感器的影响，在国控点近邻所布控的自建点上，同一时间微型空气质量检测仪所采集的数据与该国控点的数据值存在一定的差异，因此，需要利用国控点每小时的数据对国控点近邻的自建点数据进行校准。

附件1.CSV和附件2.CSV分别提供了一段时间内某个国控点每小时的数据和该国控点近邻的一个自建点数据（相应于国控点时间且间隔在5分钟内），各变量单位见附件3。请建立数学模型研究下列问题：

1. 对自建点数据与国控点数据进行探索性数据分析。
2. 对导致自建点数据与国控点数据造成差异的因素进行分析。
3. 利用国控点数据，建立数学模型对自建点数据进行校准。

2 涉及内容

2.1 涉及的技术内容

在本次实战的数据分析过程中，涉及以下技术内容：
（1）重复值的处理
（2）时间类型数据转换
（3）插值方法
（4）归一化
（5）标准化
（6）灰色关联
（7）回归模型的训练及评价
（8）各自变量重要度的评价
（9）NCA算法做特征选择
（10）模型在新数据集上预测

2.2 涉及的matlab函数

主要使用以下MATLAB功能实现：
（1）readtable
（2）unique
（3）ismissing
（4）table2array
（5）datenum
（6）interp1
（7）mapminmax
（8）mapstd
（9）fitrauto
（10）fsrnca
（11）setdiff
（12）predict

3 实战步骤

3.1 数据读取

读取数据，查看数据基本情况，

warning("off");
Data_1=readtable("附件1.csv");
Data_2=readtable("附件2.csv");
size(Data_1)
size(Data_2)

输出

3.2 去除重复值

去除时间戳重复的行，因为之后的插值和回归都属于函数问题，需保证一个时间戳只能对应一个y值。

[u_1,ia_1,ic_1]=unique(Data_1{:,end});
[u_2,ia_2,ic_2]=unique(Data_2{:,end});
u_Data_1=Data_1(ia_1,:);
u_Data_2=Data_2(ia_2,:);
size(u_Data_1)
size(u_Data_2)

输出

3.3 检查缺失值

使用函数readtable后，缺失值以NaN的形式存在，

miss_value_1=sum(sum(ismissing(u_Data_1)==true))
miss_value_2=sum(sum(ismissing(u_Data_2)==true))

输出

国建点和自建点数据均无缺失值。

3.4 异常值处理

因为我不是该项目所涉及专业领域专家，对异常值的界定标准没有经验，且我认为统计学意义上的异常并不能说明特定问题里的数据确实存在异常，又考虑到没有0值得存在，所以此处我并未对数据做处理。

3.5 初次插值尝试

%使用国建点的时间戳，对自建点数据进行插值
time1=u_Data_1{:,end};
time2=u_Data_2{:,end};
time1_num=datenum(time1);%时间类型的数据转换成数字型
time2_num=datenum(time2);%时间类型的数据转换成数字型

for i=1:6 %"两尘四气"共6列
    Data_1_spline(:,i)=interp1(time2_num,u_Data_2{:,i},time1_num,"spline");
    Data_1_linear(:,i) = interp1(time2_num,u_Data_2{:,i},time1_num,"linear");
end

输出

对比PM2.5在插值前后的表现，样条插值对较大的时间间隔和变化陡峭的点的插值较为敏感，容易出现异常的插值量，线性插值表现相对较好。

3.6 分析时间间隔

虽然线性插值表现较好，但为了进一步提升插值准确度，我打算对各个相邻时间戳的间隔进行统计分析，从而设定一个时间间隔界限值。间隔大于这个值，我将认为在这个间隔中的插值误差大，可信度低；反之，则认为可信度高，插值可取。
对自建点的时间间隔进行分析，

%分析自建点数据中两个时间点的间隔
jiange_time2_num=[];
for i=2:length(time2_num)
    jiange_time2_num(i-1)=time2_num(i)-time2_num(i-1);
end

jiange_time2_num_t=tabulate(jiange_time2_num);
jiange=jiange_time2_num_t(:,1);
bei=[];
for i=2:length(jiange)
    bei(i-1)=jiange(i)/jiange(1);%求出各个时间间隔相对于1min的倍数，也即多少分钟
end

各时间间隔频次表如下，间隔最小、频数最大的为1分钟。

3.7 有效插值时间间隔确定

接上文，为保证后续的插值的有效性，我只对时间间隔小于某个范围内的国建点时间戳进行插值，我在这里确定间隔30min为可接受的有效插值间隔。计算这时有效插值间隔数量占总时间间隔数量的比值，

sum_t_57=sum(jiange_time2_num_t(1:57,3))

得到

发现有效插值间隔占总时间间隔的99.9492%，几乎包含了绝大多数时间点。所以，我计划逐个寻找每一个time1_num(k)在time2_num里的位置m~m+1，并计算time2_num(m)到time2_num(m+1)的距离，如果该距离小于30min，则留下该点time1_num(k)进行插值。

%确定time1_num里可有效插值的点time1_num的
%第一个点在time2_num第一个点之前，所以从2开始
liu=[];
for i=2:length(time1_num)
    ind_1=find(time2_num<time1_num(i));
    ind_2=ind_1(end);
    if (time2_num(ind_2+1)-time2_num(ind_2)) <0.0208
        liu(i)=1;
    end  
end
liu(1)=1; 
sum(liu==1)/length(liu)

经过筛选，计算得出time1_num保留的时间戳数量的占time1_num总数的96.3%，绝大多数点被选中。

3.8 重新进行插值

使用筛选后的time1_num对两尘四气和其它五个影响因素都分别重新进行插值，并画图比对。由于spline方法容易出现异常点，我这里只使用linear方法。

time1_num_new=time1_num(liu==1);
Data_1_spline=[];
Data_1_linear=[];
for i=1:6 %"两尘四气"共6列
    Data_1_linear(:,i) = interp1(time2_num,u_Data_2{:,i},time1_num_new,"linear");
end

j=1;
Data_1_linear_factor=[];

for i=7:11 %影响因素共5列
    Data_1_linear_factor(:,j) = interp1(time2_num,u_Data_2{:,i},time1_num_new,"linear");
    j=j+1;
end
Data_1_linear_factor(1,:) =Data_1_linear_factor(2,:);
figure,
plot(time2_num,u_Data_2{:,4},'-ro',time1_num_new,Data_1_linear(:,4),'bo');
figure,plot(time2_num,u_Data_2{:,4},'-ro');
figure,plot(time1_num_new,Data_1_linear(:,4),'bo');

对比图，

可见由于我们设定了最小插值时间间隔，插值函数未对间隔大于30min的时间戳进行插值，尽量保证了插值的准确度，避免因为插值产生的误差在后续求关联度和进行回归时被放大。

3.9 灰色关联

先进行灰色关联前数据预处理。网上关于这个预处理，大家用的方法各有不同。有归到同一起点的，有归到平均值的（除以平均值），有归到特定范围内的（如0~1），而且算出来的关联度都互不相同。我查阅网上大家对灰色关联的解释，觉得既然灰色关联是表征两组曲线变化趋势相似度的，那归到固定范围，应该最为合适。因为我的理解，趋势即为斜率，斜率即为比例，既然只考虑比例，那大家都归到同一固定范围不就好了？所以以下我采用归到固定范围内来预处理各列数据。如果有数学大神，请给予我理论性的指导，谢谢。

Data_1_num=table2array(u_Data_1(liu==1,1:end-1));
error_data=Data_1_linear-Data_1_num;
%误差归一化
error_data_B=mapminmax(error_data',0,1);
error_data_B=error_data_B';
error_data_B(1,:)=error_data_B(2,:);
%自变量归一化
Data_X=Data_1_linear_factor;
Data_X_B=mapminmax(Data_X',0,1);
Data_X_B=Data_X_B';

灰色关联算法求关联度，

GLD=[];
for i=1:size(error_data_B,2)
    GLD(i,:)=GRA2(Data_X_B,error_data_B(:,i));  
end
GLD

关联度结果表，

所使用的灰色关联函数如下，

function out=GRA2(X,Y)
absX0_Xi=abs(X-repmat(Y,1,size(X,2)));
a=min(min(absX0_Xi));
b=max(max(absX0_Xi));
rho=0.5;
gamma=(a+rho*b)./(absX0_Xi+rho*b);
out=mean(gamma);
end

3.10 回归模型训练

先进行所需数据预处理。首先是数据标准化，当样本量足够大的时候，随机划分的训练集和测试集其实是有相同的分布的，这里在划分数据集前进行数据标准化，

[y_data,ps_1]=mapstd(error_data');%因为要向训练集范围外回归，所以未使用归一化
y_data=y_data';

[x_data,ps_2]=mapstd(Data_X');
x_data=x_data';

构建回归模型所需数据，先取误差的第一列，即将自建点和国建点PM2.5的差值作为回归目标，进行训练，

tbl_1=array2table([y_data(:,1),x_data]);
tbl_1.Properties.VariableNames=["Y","x1","x2","x3","x4","x5"];

划分数据集，

rng('default') % For reproducibility of the data partition
c = cvpartition(size(tbl_1,1),'Holdout',0.2);
trainingIdx = training(c); % Training set indices
tbl_1_Train = tbl_1(trainingIdx,:);
testIdx = test(c); % Test set indices
tbl_1_Test = tbl_1(testIdx,:);

训练回归模型，使用了fitrauto，该函数可以自动选择最优模型及对应的超参数，是我等懒人“居家旅行，摸鱼放松”的必备函数，(•‾̑⌣‾̑•)

options = struct('UseParallel',true);
Mdl = fitrauto(tbl_1_Train,'Y','OptimizeHyperparameters','all', ...
    'HyperparameterOptimizationOptions',options,...
    "Learners","tree");

这里我只选定使用树模型，便于后续对各特征打分。大家也可以尝试不同的模型，或者采用大人的方法：“我都要！”

3.11 评价模型性能

评价模型在测试集上的性能，

testMSE = loss(Mdl,tbl_1_Test,'Y');
testError = log(1 + testMSE)

模型对各自变量在回归模型中的重要性进行打分，

predictorImportance(Mdl)

使用NCA算法做特征选择，也可以得出各自变量在回归时的重要性排序，

X=tbl_1_Train{:,2:end};
y=tbl_1_Train{:,1};
N=length(y);
mdl_nca = fsrnca(X,y,'Verbose',1,'Lambda',0.5/N);
mdl_nca.FeatureWeights

可以看出，两种对特征重要度的评价方法都对第一个特征打了最低分。我觉得可以结合上述灰色关联度，对各影响因素进行综合评价。

3.12 在新数据上做预测

使用训练好的模型在新数据上做预测，这里相当于“外推”，向外拟合，

[~,ia] = setdiff(time2_num,time1_num_new);
y_base=u_Data_2{ia,1};
predict_x=u_Data_2{ia,7:11};
predict_x=mapstd('apply',predict_x',ps_2);
predict_x=predict_x';
Yfit = predict(Mdl,predict_x);
%逆标准化
Yfit = mapstd('reverse',Yfit',ps_1);
Yfit =Yfit';
y_new=y_base-Yfit;

画图比对，

figure,
plot(time1_num,u_Data_1{:,1},'ro','LineWidth',2),hold on;
plot(time2_num(ia),y_new(:,1),'bo');hold off;

figure,
plot(time1_num,u_Data_1{:,1},'ro','LineWidth',2),hold on;
plot(time2_num,u_Data_2{:,1},'go');hold off;

可见，未修正前的自建点数据（图二），整体要比国建点偏上，而修正后的自建点数据（图一），整体将国建点“包裹”的较好，更加合理。当然，使用回归的方法修正后，自建点出现了一些负值，将这些负值置0即可。
大家可以尝试对剩下的“一尘四气”进行回归，并比对效果，此处不赘述了。

来源：木星流火