1 題目簡介

題目來源于2019 高教社杯全國大學生數學建模競賽D題——空氣質量數據的校準。
空氣污染對生態環境和人類健康危害巨大，通過對“兩塵四氣”（PM2.5、PM10、CO、NO2、SO2、O3）濃度的實時監測可以及時掌握空氣質量，對污染源采取相應措施。雖然國家監測控制站點（國控點）對“兩塵四氣”有監測數據，且較為準確，但因為國控點的布控較少，數據發布時間滯后較長且花費較大，無法給出實時空氣質量的監測和預報。某公司自主研發的微型空氣質量檢測儀（如圖所示）花費小，可對某一地區空氣質量進行實時網格化監控，并同時監測溫度、濕度、風速、氣壓、降水等氣象參數。

由于所使用的電化學氣體傳感器在長時間使用后會產生一定的零點漂移和量程漂移，非常規氣態污染物（氣）濃度變化對傳感器存在交叉干擾，以及天氣因素對傳感器的影響，在國控點近鄰所布控的自建點上，同一時間微型空氣質量檢測儀所采集的數據與該國控點的數據值存在一定的差異，因此，需要利用國控點每小時的數據對國控點近鄰的自建點數據進行校準。

附件1.CSV和附件2.CSV分別提供了一段時間內某個國控點每小時的數據和該國控點近鄰的一個自建點數據（相應于國控點時間且間隔在5分鐘內），各變量單位見附件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); endjiange_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.0208liu(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"); endj=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即可。
大家可以嘗試對剩下的“一塵四氣”進行回歸，并比對效果，此處不贅述了。

總結

以上是生活随笔為你收集整理的实战1 - 空气质量数据的校准的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。