2019年国赛高教杯数学建模
D题 空气质量数据的校准
原题再现
空气污染对生态环境和人类健康危害巨大，通过对“两尘四气”（PM2.5、PM10、CO、NO2、SO2、O3）浓度的实时监测可以及时掌握空气质量，对污染源采取相应措施。虽然国家监测控制站点（国控点）对“两尘四气”有监测数据，且较为准确，但因为国控点的布控较少，数据发布时间滞后较长且花费较大，无法给出实时空气质量的监测和预报。某公司自主研发的微型空气质量检测仪（如图所示）花费小，可对某一地区空气质量进行实时网格化监控，并同时监测温度、湿度、风速、气压、降水等气象参数。
由于所使用的电化学气体传感器在长时间使用后会产生一定的零点漂移和量程漂移，非常规气态污染物（气）浓度变化对传感器存在交叉干扰，以及天气因素对传感器的影响，在国控点近邻所布控的自建点上，同一时间微型空气质量检测仪所采集的数据与该国控点的数据值存在一定的差异，因此，需要利用国控点每小时的数据对国控点近邻的自建点数据进行校准。
附件 1.CSV 和附件 2.CSV 分别提供了一段时间内某个国控点每小时的数据和该国控点近邻的一个自建点数据（相应于国控点时间且间隔在 5 分钟内），各变量单位见附件 3。请建立数学模型研究下列问题：
1. 对自建点数据与国控点数据进行探索性数据分析。
2. 对导致自建点数据与国控点数据造成差异的因素进行分析。
3. 利用国控点数据，建立数学模型对自建点数据进行校准。
整体求解过程概述(摘要)
掌握空气质量并采取相应措施，对维持生态环境和推动人类健康具有重大的现实意义。本文以空气质量的监测数据为研究对象，在合理假设的基础上，首先，对国控点和自建点的监测数据进行了探索性分析；其次，找到了一种基于时间对齐的差异值计算方法，在此基础上分别建立并求解了基于关联度分析和基于多元线性回归的差异值影响模型，分析了导致自建点数据和国控点数据产生差异的因素；接着，基于分段线性插值的自建点数据校准模型解答了问题三；最后，对模型优缺点进行了分析。主要结论如下：
针对问题一，从数据质量检查、数据的统计量、污染物浓度直方图、污染物浓度和气象参数的相关系数、“两尘四气”浓度的 24 小时波动趋势等五个方面对自建点数据和国控点数据开展探索性分析。结论表明：国控点数据中存在多处“小时数的不连续点”，自建点数据也存在一些不连续点；该区域空气污染主要由颗粒物造成，PM2.5、PM10达到二级浓度限值，“四气”均达到一级浓度限值；从直方图看，国控点和自建点数据在分布上存在明显差异； “两尘四气”浓度和气象参数等 11 个因素之间存在显著的两两相关关系；“两尘四气”浓度在一天内达到峰值与谷值的时间段不同，国控点和自建点数据污染物浓度的 24 小时变动趋势不一致。
针对问题二，基于时间对齐算法实现国控点数据和自建点数据的可比较性；建立并求解了基于关联度分析和多元线性回归的差异值影响模型，均循环 6 次，分别考察了 11个影响因素对“两尘四气”浓度差异值的影响。时间对齐后，附件 1 的有效记录为 4137行，附件 2 的有效记录为 197341 行，进而将附件 2 压缩为 4137 行，定量计算出两者的差异值。关联度分析结论表明：影响“两尘四气”浓度差异值的最主要因素不同，对 PM2.5、PM10、CO、NO2、SO2、O3 监测浓度差异影响最大的分别为第 11 个（湿度）、第 11 个（湿度）、第 4 个（NO2）、第 5 个（SO2）、第 5 个（SO2）、第 6 个（O3）影响因素；污染物浓度变化的确会对传感器存在交叉干扰；影响因素的重要程度并非一成不变，针对不同污染物，影响因素的主次关系变化较大。通过多元回归分析发现：P 值表明所有 6个模型在整体上是显著的，R2 值表明模型 2、模型 3 和模型 6 的拟合度较高；零点漂移和量程漂移可视作回归方程的常数项变化；传感器的长时间使用、气体污染物浓度变化以及气象因素均可导致自建点数据和国控点数据之间的差异，影响程度由回归系数决定；并非只有气体污染物对数据差异产生影响，颗粒物的影响效应同样存在。
针对问题三，基于分段线性插值方法对自建点数据进行校准，将自建点数据校准问题，转化为一个过已知有限个数据点（国控点监测数据）求近似函数的问题，从而将附件 1 中的 4137 条国控点监测数据，扩展成为 197341 条自建点的校准数据。结论表明：校准后的自建点数据与国控点数据吻合度较好；自建点与国控点的残差为 56717、校准后的数据与国控点的残差为 1649，仅为原残差的 2.9%。文中所建立的模型简便易行，便于推广，可利用国控点的数据对近邻自建点的数据进行校准。
模型假设：
为简化问题，做出如下合理假设：
（1）同一时间，自建点与国控点所处的客观环境是一致的，即污染物浓度和天气因素的实际值是完全相同的。
（2）不考虑监测数据的存储、传输、分发等环节的随机噪声所导致的数据误差。
（3）假设数据总体是符合线性正态误差模型，即误差项 ε ~ N(0,σ^2)
问题分析：
考虑到数据量较大，首先需要对附件 1.CSV（国控点监测数据），附件 2.CSV（自建点监测数据，含气象参数）进行数据预处理，将后续程序中需要用到的数据另存为后缀为.mat 的文件。
探索性数据分析（Exploratory Data Analysis，简称 EDA），是指对已有的数据（题目中对应为附件1和附件2中的原始数据）在尽量少的先验假定下进行探索，通过作图、制表、方程拟合、计算特征量等手段探索数据的结构和规律的一种数据分析方法[2]。针对问题一，首先需要对原始数据进行数据质量检查，检查内容主要包括数据一致性，处理无效值和缺失值；在此基础上，通过最小值、最大值、均值、中位数、标准差、偏度以及峰度等统计量反映国控点和自建点监测数据的数量特征；进一步地，还可绘制出变量的直方图、分析变量之间的相关关系，通过可视化的方式直观考察国控点和自建点在监测值上的差异性。
影响差异值的三大来源分别是：传感器的零点漂移和量程漂移、污染物的交叉干扰、天气因素对传感器的影响。针对问题二，主要包括两个子问题：1）如何定量地刻画国控点和自建点监测数据的差异值；2）如何建立合理模型，对引致差异值的影响因素进行分析。分析要点应该包括：究竟是哪些影响因素导致了差异值的产生？能否通过分析关联度，测算出影响因素的排序？是否可以通过多元回归，考察各影响因素对于差异值的显著程度、作用方向及影响程度？文中考虑采用灰度系统理论中的关联度分析方法和多元线性回归方法。
针对问题三，已知国控点的监测数据准确，但布控较少、发布时间滞后；自建点监测数据更新快，但误差较大。如何利用国控点数据对自建点数据进行校准这一问题，等同于求一个过已知有限个数据点（国控点监测数据）的近似函数，进而产生出与自建点监测数据同步更新的数据，并据此对自建点数据进行校准。校准方法的有效性需要进行客观、准确的评价，可考虑通过可视化方式、定量计算两种方式进行评价。文中考虑采用线性插值的方法进行数据校准。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：(代码和文档not free)
The actual procedure is shown in the screenshot
function data = getTime(str,flag)
% 函数功能：根据输入的时间字符串，flag=2,输出相应的年、月、日
% flag 为其它数值时，输出相应的年、月、日、时、分、秒
if nargin ==1
flag = 2;
end
tmpdate = datevec(str);
if flag == 1
data = tmpdate(1:3,:);
else
data = tmpdate;
end
% 对数据进行预处理
% Tab1 预处理：标题行 tab1Title、国控点时间 tab1Time、监测值 tab1MonitorValue
[~,~,tab1Alldata]= xlsread('附件 1.csv'); % 读取“附件 1.csv”
Ltab1 = length(tab1Alldata); % 取行数
% tab1Title = tab1Alldata(1 ,:); % tab1 行名称
tab1Alldata = tab1Alldata(2:Ltab1 ,:); % 删除标题行
tab1Time = getTime(tab1Alldata(:,7)); % 国控点时间：年月日时分秒
tab1Time = tab1Time(:,1:4); % 仅保留年、月、日、时
tab1MonitorValue = cell2mat(tab1Alldata(:,1:6)); % tab1 传感器记录值
% save tab1Title.mat tab1Title
save tab1Time.mat tab1Time
save tab1MonitorValue.mat tab1MonitorValue
% Tab2 预处理：标题行 tab2Title、自建点时间 tab2Time、监测值 tab2MonitorValue
[~,~,tab2Alldata]= xlsread('附件 2.csv'); % 读取“附件 2.csv”
Ltab2 = length(tab2Alldata); % 取行数
% tab2Title = tab2Alldata(1 ,:); % tab2 行名称
tab2Alldata = tab2Alldata(2:Ltab2 ,:); % 删除标题行
tab2Time = getTime(tab2Alldata(:,12)); % 自建点时间：年月日时分秒
tab2Time = tab2Time(:,1:5); % 仅保留年、月、日、时、分
tab2MonitorValue = cell2mat(tab2Alldata(:,1:11)); % tab2 传感器记录值
% save tab2Title.mat tab2Title
save tab2Time.mat tab2Time
save tab2MonitorValue.mat tab2MonitorValue
% 解决问题 1，对自建点数据和国控点数据进行探索性数据分析
%% 装载数据
% 国控点数据 附件 1 标题：PM2.5、PM10、CO、NO2、SO2、O3
load tab1Time.mat % 附件 1 国控点记录时间：年月日时 （4200 * 4）
load tab1MonitorValue.mat % 附件 1 监测值：4200 * 6
% 自建点数据
% 附件 2 标题：PM2.5、PM10、CO、NO2、SO2、O3、风速、压强、降水量、温度、湿度
load tab2Time.mat % 附件 2 自建点记录时间：年月日时分 （234717 * 5）
load tab2MonitorValue.mat % 附件 2 监测值：（234717 * 11）
%% 缺失值分析
% ii = 1:4200;plot(ii,tab1Time(:,4));ii = 1:1000; plot(ii,tab1Time(ii,4));
ii = 300:600;plot(ii,tab1Time(ii,4),'-b','LineWidth',2);grid on
xlabel('附件 1 的数据记录号','FontSize',12);
ylabel('记录时间的小时数','FontSize',12);
set(gca,'ytick',0:2:23)
set(gca,'FontName','Times New Roman','FontSize',14) %设置坐标轴刻度字体名称、大小
%% 描述性统计：最小值、最大值、均值、中位数、标准差、偏度、峰度
% 国控点数据
tab1DataDisplay = [];
for i = 1:6
tmpDataDisplay = explorDataAnalysis(tab1MonitorValue(:,i));
tab1DataDisplay = [tab1DataDisplay tmpDataDisplay];
end
save tab1DataDisplay.mat tab1DataDisplay
% 自建点数据
tab2DataDisplay = [];
for i = 1:11
tmpDataDisplay = explorDataAnalysis(tab2MonitorValue(:,i));
tab2DataDisplay = [tab2DataDisplay tmpDataDisplay];
end
save tab2DataDisplay.mat tab2DataDisplay
%% 相关系数：结果与 corrcoef 相同
C1 = corr(tab1MonitorValue,'type','Pearson'); % 国控点数据
C2 = corr(tab2MonitorValue,'type','Pearson'); % 自建点数据
%% 国控点数据波形图
figure
tab1Date = 1:length(tab1Time);
set(gca,'FontName','Times New Roman','FontSize',12) %设置坐标轴刻度字体名称、大小
subplot(231);plot(tab1Date,tab1MonitorValue(:,1));title('国控点 PM2.5','FontSize',12);xlim([0 4300]);
subplot(232);plot(tab1Date,tab1MonitorValue(:,2));title('国控点 PM10','FontSize',12);xlim([0 4300]);
subplot(233);plot(tab1Date,tab1MonitorValue(:,3));title('国控点 CO','FontSize',12);xlim([0 4300]);
subplot(234);plot(tab1Date,tab1MonitorValue(:,4));title('国控点 NO2','FontSize',12);xlim([0 4300]);
subplot(235);plot(tab1Date,tab1MonitorValue(:,5));title('国控点 SO2','FontSize',12);xlim([0 4300]);
subplot(236);plot(tab1Date,tab1MonitorValue(:,6));title('国控点 O3','FontSize',12);xlim([0 4300]);
%% 对 24 小时浓度均值的波动情况进行统计分析
% 国控点数据
hourWave = zeros(24,6);
for i = 0:23
tmpIndex = tab1Time(:,4)==i;
tmpValue = tab1MonitorValue(tmpIndex,:);
hourWave(i+1,:) = sum(tmpValue,1)/length(tmpValue);
end
% 自建点数据
hourWave2 = zeros(24,6);
for i = 0:23
tmpIndex = tab2Time(:,4)==i;
tmpValue = tab2MonitorValue(tmpIndex,1:6);
hourWave2(i+1,:) = sum(tmpValue,1)/length(tmpValue);
end
figure
tab1Date = 0:23;
set(gca,'FontName','Times New Roman','FontSize',12,'LineWidth',3)%设置坐标轴刻度字体、大小、线宽
subplot(231);plot(tab1Date,hourWave(:,1),'--bs','MarkerFaceColor','g');
hold on;plot(tab1Date,hourWave2(:,1),'-ro');
title('PM2.5 均值 24 小时波动','FontSize',12);xlim([0 24]);
subplot(232);plot(tab1Date,hourWave(:,2),'-bs','MarkerFaceColor','g');
hold on;plot(tab1Date,hourWave2(:,2),'-ro');
title('PM10 均值 24 小时波动','FontSize',12);xlim([0 24]);
subplot(233);plot(tab1Date,hourWave(:,3),'-bs','MarkerFaceColor','g');
hold on;plot(tab1Date,hourWave2(:,3),'-ro');
title('CO 均值 24 小时波动','FontSize',12);xlim([0 24]);
subplot(234);plot(tab1Date,hourWave(:,4),'-bs','MarkerFaceColor','g');
hold on;plot(tab1Date,hourWave2(:,4),'-ro');
title('NO2 均值 24 小时波动','FontSize',12);xlim([0 24]);
subplot(235);plot(tab1Date,hourWave(:,5),'-bs','MarkerFaceColor','g');
hold on;plot(tab1Date,hourWave2(:,5),'-ro');
title('SO2 均值 24 小时波动','FontSize',12);xlim([0 24]);
subplot(236);plot(tab1Date,hourWave(:,6),'-bs','MarkerFaceColor','g');
hold on;plot(tab1Date,hourWave2(:,6),'-ro');
title('O3 均值 24 小时波动','FontSize',12);xlim([0 24]);
legend('国控点','自建点','FontSize',12)
%% tab1 和 tab2 直方图比较
figure
binNum = 10;
set(gca,'FontName','Times New Roman','FontSize',12) %设置坐标轴刻度字体名称、大小
subplot(341);hist(tab1MonitorValue(:,1),binNum);title('国控点 PM2.5 直方图','FontSize',12);
subplot(342);hist(tab2MonitorValue(:,1),binNum);title('自建点 PM2.5 直方图','FontSize',12);
subplot(343);hist(tab1MonitorValue(:,2),binNum);title('国控点 PM10 直方图','FontSize',12);
subplot(344);hist(tab2MonitorValue(:,2),binNum);title('自建点 PM10 直方图','FontSize',12);
subplot(345);hist(tab1MonitorValue(:,3),binNum);title('国控点 CO 直方图','FontSize',12);
subplot(346);hist(tab2MonitorValue(:,3),binNum);title('自建点 CO 直方图','FontSize',12);
subplot(347);hist(tab1MonitorValue(:,4),binNum);title('国控点 NO2 直方图','FontSize',12);
subplot(348);hist(tab2MonitorValue(:,4),binNum);title('自建点 NO2 直方图','FontSize',12);
subplot(349);hist(tab1MonitorValue(:,5),binNum);title('国控点 SO2 直方图','FontSize',12);
subplot(3,4,10);hist(tab2MonitorValue(:,5),binNum);title('自建点 SO2 直方图','FontSize',12);
subplot(3,4,11);hist(tab1MonitorValue(:,6),binNum);title('国控点 O3 直方图','FontSize',12);
subplot(3,4,12);hist(tab2MonitorValue(:,6),binNum);title('自建点 O3 直方图','FontSize',12);
642
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
