模式+机器学习如何用于气候变化预估：一个海温订正与预测案例EEMD-BPNN

本文节选自和鲸社区《Python气象海洋地理数据分析Workshop》第一期教案，报名即可获取完整教案及更多精彩内容；教案使用和鲸ModelWhale工具撰写，支持一键跑通，点下方报名 👇

导师：宋振亚 自然资源部第一海洋研究所研究员、博导

助教：匡志远 宋振亚老师团队成员，教案主要撰写者

 · 简介 · 

本期，我们将基于和鲸 ModelWhale 平台，手把手教大家动手开展气候模式预估全球平均海表面温度订正，并对模式结果进行后处理。我们将采用数据分解和机器学习结合的方法来订正气候模式模拟和预估的全球平均海表温度偏差。

 · 目录 · 

SST数据预处理
1）读取数据
2）空间插值
3）计算全球平均SST
4）可视化

SST数据的EEMD分解
1）EEMD方法的工作原理

EEMD后处理
1）组合时间序列以满足物理约束

SST订正及预估
1）海表面温度偏差订正模型
2）数据处理
3）评价订正结果

作业
1）结合项目教学内容，进行模型调优
2）尝试用同样的方法，订正另一个数据源数据

 · SST数据预处理 · 

数据简介

对SST数据预处理之前，我们需要先对数据有一个大致的了解。本章节我们选用的数据是自然资源部第一海洋研究所自主发展的第二代气候模式FIO-ESM v2.0模拟和预估的全球SST以及观测值，模式和观测的数据格式都是nc格式。

选用的数据为：
1） 观测数据：选用的是第5代扩展重构SST数据ERSST v5 ，空间覆盖全球海洋，分辨率为2°×2°，时间范围从1854年1月到2019年12月，时间间隔是1个月，共166年的数据。
2） 模式数据：包含历史时期模拟数据和未来情景预估数据。FIO-ESM v2.0历史时期数据，空间覆盖全球海洋，高纬度地区1.1°，赤道地区加密为0.3°—0.5°，时间范围从1850年1月到2014年12月， 时间间隔是1个月，共165年的数据; 3种未来不同温室气体排放程度预估的数据，空间覆盖全球海洋，分辨率和历史数据相同，时间范围从2015年1月到2100年12月，时间间隔是1个月，共86年的数据。

读取数据

略

观测数据插值

观测数据和模式数据的网格分辨率不相同，特别是模式数据还是不均匀的网格，所以不利于我们后面对其计算全球平均，所以我们为了统一，先将观测和模式数据进行空间插值到相同网格下，方便后面的计算。

这里，我们使用griddata函数对数据进行插值，具体插值方法介绍可见python的griddata插值。

## 插值前的经纬度, 维度是(180*89,2)
LatLon_ERSST_Before = np.hstack( ( Lat_ERSST_mgrid.reshape(-1,1), Lon_ERSST_mgrid.reshape(-1,1) ) ) 
## 插值前的SST，维度是(180*89,1)
SST_ERSST_Before = SST_ERSST.reshape( (SST_ERSST.shape[0],-1,1) )
## 插值前的经纬度和SST有了，我们还需要给出插值后的经纬度信息
Lat_ERSST_After_1D = np.arange(Max_Lat_ERSST,Min_Lat_ERSST-1,-Step_LatLon_After) 
Lon_ERSST_After_1D = np.arange(Min_Lon_ERSST,Max_Lon_ERSST+1,Step_LatLon_After)
Lat_ERSST_After_2D, Lon_ERSST_After_2D = np.meshgrid(Lat_ERSST_After_1D,Lon_ERSST_After_1D)
Lat_ERSST_After_2D = Lat_ERSST_After_2D.T
Lon_ERSST_After_2D = Lon_ERSST_After_2D.T
## 使用griddata函数对SST插值
SST_ERSST_After = np.zeros( [SST_ERSST.shape[0], Lat_ERSST_After_2D.shape[0],  Lat_ERSST_After_2D.shape[1] ])
## 对每一个月的全球观测资料插值，并⏲计时
start_time = time.time()
for t in range(SST_ERSST.shape[0]):
    # griddata函数插值，需要输入插值前的经纬度信息，插值前的变量，插值后的经纬度，插值方法，这里插值方法选择'nearest'方法。
    SST_ERSST_After_Now = griddata(LatLon_ERSST_Before, SST_ERSST_Before[t,:,:], (Lat_ERSST_After_2D, Lon_ERSST_After_2D), method='nearest').squeeze()
    SST_ERSST_After[t, :, :] = SST_ERSST_After_Now
    if t % 100 == 0


    
: # 每100个数打印显示一下
        print(t)
end_time = time.time()
print(str(end_time - start_time) + "s") # 打印看一下循环计算的用时

> 89.94760847091675s

模式的历史数据插值、模式的未来预估数据插值，操作与上面类似。

计算全球平均SST

这里计算全球平均我们要注意不同纬度下的网格面积是不一样的，所以需要采用加权平均，不能用算术平均，否则会带来较大误差。

## 定义求区域平均函数
def areamean2D(X,lat,lon):
    """
    本函数计算经纬度网格下某变量区域平均值。
    
    输入：
        X: 待求区域平均的变量，如SST;
        lat: 纬度；
        lon: 经度。
    输出：
        X_AreaMean: 输入变量X的区域平均值
 
    """
    import numpy as np
 
    Len_lat = lat.shape[0]
    Len_lon = lon.shape[0]
    lat_2D = lat.repeat(Len_lon).reshape((Len_lat,Len_lon))
    Weight = np.cos(np.pi/180*lat_2D)  # 关键是要计算不同纬度对应的权重
    Weight_Ocean = Weight[ np.where( abs(X) 50 ) ]  # 设置一个较大的较为合理的值，小于这个值才是正常的SST值
    X_Ocean = X[ np.where( abs(X) 50 ) ]
    X_AreaMean = sum( X_Ocean * Weight_Ocean ) / sum( Weight_Ocean )
 
    return X_AreaMean  # 返回计算结果

计算观测的全球平均SST

SST_ERSST_AreaMean = np.zeros(SST_ERSST_After.shape[0]) # 预定义一个空数组，以储存观测的全球平均SST。
## 开始对观测SST逐月计算全球平均值
start_time = time.time()
for t in range(SST_ERSST_After.shape[0]):
    # 调用函数areamean2D
    SST_ERSST_AreaMean[t] = areamean2D(SST_ERSST_After[t,:,:], Lat_ERSST_After_1D, Lon_ERSST_After_1D)
    if t % 100 == 0: # 每100个数打印显示一下
        print(t)
end_time = time.time()
print(str(end_time - start_time) + "s")

计算模式历史的全球平均SST、模式未来预估的全球平均SST，操作与上面类似

绘图

历史时期全球平均SST的观测、模式和模式偏差的时间序列（模式减观测）

上图自上到下依次是：a）观测SST，b）模式SST，c）模式模拟的SST偏差（模式减去观测）在185401-201412的全球平均时间序列，可以看出观测和模式的全球平均SST随着时间均呈现增长的趋势，但是模式相对观测还存在正负1度之间的偏差，且偏差时间序列有着非线性非平稳的特征。

模式预估全球平均SST的3种未来情景试验的时间序列

上图自上到下依次是：a）低排放预估，b）中等排放预估，c）高排放预估在201501-210012的全球平均时间序列，三种排放预估结果都显示未来全球平均SST将随时间逐渐升高，排放越高，上升得越快，到本世纪末温度越高。

 · SST数据的EEMD分解 · 

得到全球平均的SST时间序列之后，我们先不直接拿来做订正，因为原始时间序列长度是百年尺度的非线性非平稳序列，本身会包含多种时间尺度（如季节、年、十年、几十年等）的信息，直接订正效果未必好。如果将原始时间序列分解成多种时间尺度的序列，就相当于得到多种更为丰富的时间序列信息，同时有了不同时间尺度的物理约束。

基于上述的想法，我们采用EEMD（Ensemble Empirical Mode Decomposition）方法对原始时间序列来做，EEMD能将原始数据分解成多个从高频到低频的序列，也称本征模函数IMF（Intrinsic Mode Function），分解后的IMF加和也能以可忽略的误差范围接近原始序列。我们可以按照不同时间尺度对IMF组合，将各个组合时间序列输入到机器学习模型里做订正，最后将订正后的组合加起来得到最终的订正时间序列，这样比原始序列直接订正会更精确，且有较清楚的物理约束。

所以，这一章，我们的目的很明确，就是定义EEMD函数，然后调用函数去分解观测和模式数据。

定义EEMD函数

这里我们需要定义两个用于EEMD数据分解的函数，一个是求极值的函数”extrema“，另一个是主函数—“eemd”函数。大家在调用函数的时候只需要输入两个参数和待分解的变量，里面的细节感兴趣可以自行学习。

EEMD分解观测和模式SST 

eemd分解除了输入原始时间序列，还需要输入两个参数，一个是加入噪声的标准差与输入数据标准差的比值，也就是Nstd，一般在0-0.2之间取值。另一个是 EEMD使用的集合数，也就是NE，常取50，100，200等值。这里，对观测和模式数据我们均设置相同的值，即Nstd取0.2，NE取100。

% 设置两个参数
Nstd = 0.2;
NE = 100;
% 开始对观测数据分解，并计时
tic
EEMD_ERSST = eemd( SST_ERSST_AreaMean, Nstd, NE ); % 得到第一列是原始序列，第二列到最后一列是分解的IMF分量
toc  
IMFs_ERSST = EEMD_ERSST(:, 2:end)';
NIMFs_ERSST= size(IMFs_ERSST,1);
Recons_IMFs_ERSST = sum(IMFs_ERSST, 1);
Delta_IMFs_ERSST = SST_ERSST_AreaMean - Recons_IMFs_ERSST;
print((Delta_IMFs_ERSST0, max(Delta_IMFs_ERSST)) % 可以看看残差的最大值最小值
% 画图，可以先不画，后面我们把数据保存好了再下一章python3环境下画。
% t = 1:Num_Sample_SST;
% hf = figure;
% for n = 1:NIMFs_ERSST
%   subplot(NIMFs_ERSST, 1, n);
%    plot(t, IMFs_ERSST(n,:) );
%    ylabel( ['IMF', num2str(n) ] );
% end 
% xlabel( 'Time' );

分解模式的历史SST、分解模式的未来预估SST，操作与上面类似

 · EEMD后处理 · 

画图：画出观测和模式历史和未来的各个IMF

画观测SST的IMF

上图是观测SST在192901-201412的10个IMF时间序列，自上到下依次是IMF1到IMF10。其中IMF1到IMF9频率依次降低，IMF10是趋势项，可以看出近百年来，我们生存的地球平均SST呈现逐年升温的趋势。

画模式的历史SST的IMF

上图是模式SST在192901-201412的10个IMF时间序列，自上到下依次是IMF1到IMF10。其中IMF1到IMF9频率依次降低，IMF10是趋势项，可以看出近百年来，模式也能模拟出全球平均SST呈现逐年升温的趋势。

画模式的未来预估SST的IMF

上图是模式在未来低排放情景下预估SST在201501-210012的10个IMF时间序列，自上到下依次是IMF1到IMF10。其中IMF1到IMF9频率依次降低，IMF10是趋势项，预估结果显示，未来近百年全球平均SST呈现继续升温趋势。

SSP245

上图是模式在未来中等排放情景下预估SST在201501-210012的10个IMF时间序列，自上到下依次是IMF1到IMF10。其中IMF1到IMF9频率依次降低，IMF10是趋势项，预估结果显示，未来近百年全球平均SST呈现继续升温趋势，升温幅度比低排放更大。

SSP585

上图是模式在未来高排放情景下预估SST在201501-210012的10个IMF时间序列，自上到下依次是IMF1到IMF10。其中IMF1到IMF9频率依次降低，IMF10是趋势项，预估结果显示，未来近百年全球平均SST呈现继续升温趋势，升温幅度在三种排放情景中最大，本世纪末较2015年升高近4°C。

计算EEMD分解数据平均周期

计算各个IMF的平均周期，大家要特别注意最后两个IMF，计算公式可能带来较大误差

# ERSST
Num_CrossZero_ERSST = np.zeros(NIMFs_ERSST - 1)
MeanPeriod_IMFs_ERSST = np.zeros(NIMFs_ERSST - 1)
MeanYearPeriod_IMFs_ERSST = np.zeros(NIMFs_ERSST - 1)
for n in range(NIMFs_ERSST - 1):
    k = 0
    for i in range(len(Time_SST) - 1):
        if


    
 IMFs_ERSST[n,i] * IMFs_ERSST[n,i + 1] 0:  # 数跨零点的个数
            k = k + 1
    Num_CrossZero_ERSST[n] = k
    MeanPeriod_IMFs_ERSST[n] = len(Time_SST) * 2 / k  # 根据跨零点的个数计算平均周期，单位是月
    MeanYearPeriod_IMFs_ERSST[n] = len(Time_SST) * 2 / k / 12  # 单位转换成年
print(MeanYearPeriod_IMFs_ERSST)

# Historical
NIMFs_Historical = IMFs_Historical.shape[0]
Num_CrossZero_Historical = np.zeros(NIMFs_Historical - 1)
MeanPeriod_IMFs_Historical = np.zeros(NIMFs_Historical - 1)
MeanYearPeriod_IMFs_Historical = np.zeros(NIMFs_Historical - 1)
for n in range(NIMFs_Historical - 1):
    k = 0
    for i in range(len(Time_SST)-1):
        if IMFs_Historical[n,i] * IMFs_Historical[n,i + 1] 0:
            k = k + 1
    Num_CrossZero_Historical[n] = k
    MeanPeriod_IMFs_Historical[n] = len(Time_SST) * 2 / k
    MeanYearPeriod_IMFs_Historical[n] = len(Time_SST) * 2 / k / 12
print(MeanYearPeriod_IMFs_Historical)

# SSP126/SSP245/SSP585
NIMFs_SSP = IMFs_SSP126.shape[0]
Num_CrossZero_SSP126 = np.zeros(NIMFs_SSP - 1)
Num_CrossZero_SSP245 = np.zeros(NIMFs_SSP - 1)
Num_CrossZero_SSP585 = np.zeros(NIMFs_SSP - 1)
MeanPeriod_IMFs_SSP126 = np.zeros(NIMFs_SSP - 1)
MeanPeriod_IMFs_SSP245 = np.zeros(NIMFs_SSP - 1)
MeanPeriod_IMFs_SSP585 = np.zeros(NIMFs_SSP - 1)
MeanYearPeriod_IMFs_SSP126 = np.zeros(NIMFs_SSP - 1)
MeanYearPeriod_IMFs_SSP245 = np.zeros(NIMFs_SSP - 1)
MeanYearPeriod_IMFs_SSP585 = np.zeros(NIMFs_SSP - 1)
for n in range(NIMFs_SSP - 1):
    # SSP126
    k = 0
    for i in range(len(Time_SST)-1):
        if IMFs_SSP126[n,i] * IMFs_SSP126[n,i + 1] 0:
            k = k + 1
    Num_CrossZero_SSP126[n] = k
    MeanPeriod_IMFs_SSP126[n] = len(Time_SST) * 2 / k
    MeanYearPeriod_IMFs_SSP126[n] = len(Time_SST) * 2 / k / 12

    # SSP245
    k = 0
    for i in range(len(Time_SST) - 1):
        if IMFs_SSP245[n, i] * IMFs_SSP245[n, i + 1] 0:
            k = k + 1
    Num_CrossZero_SSP245[n] = k
    MeanPeriod_IMFs_SSP245[n] = len(Time_SST) * 2 / k
    MeanYearPeriod_IMFs_SSP245[n] = len(Time_SST) * 2 / k / 12

    # SSP585
    k = 0
    for i in range(len(Time_SST) - 1):
        if IMFs_SSP585[n, i] * IMFs_SSP585[n, i + 1] 0:
            k = k + 1
    Num_CrossZero_SSP585[n] = k
    MeanPeriod_IMFs_SSP585[n] = len(Time_SST) * 2 / k
    MeanYearPeriod_IMFs_SSP585[n] = len(Time_SST) * 2 / k / 12

print(MeanYearPeriod_IMFs_SSP126)
print(MeanYearPeriod_IMFs_SSP245)
print(MeanYearPeriod_IMFs_SSP585)

> [ 0.49283668  0.89583333  1.15436242  3.51020408  7.16666667 13.23076923 43.         43.         86.        ] [ 0.48450704  1.00584795  1.49565217  4.41025641  7.47826087 19.11111111 43.         86.         86.        ] [ 0.48450704  1.00584795  1.4214876   4.77777778  6.88       15.63636364 43.         86.         86.        ] [ 0.45382586  0.86432161  1.00584795  3.51020408  4.91428571 11.46666667 24.57142857 43.         86.        ] [ 0.37310195  0.55844156  1.00584795  2.60606061  4.64864865  9.05263158 43.         86.         86.        ]

根据平均周期组合IMF 我们依据的时间尺度有以下6个:

1. 季节：3个月~12个月
2. 年：12个月左右
3. 年际：1年~10年
4. 十年：10年左右
5. 年代际：大于10年
6. 年代际：>80年

 · SST订正及预估 · 

BPNN机器学习模型简介

BPNN（Back Propagation Neural Network）一般包括输入层、隐含层和输出层。一般地，输入层和输出层只有一层，隐含层则至少有一层。网络的每一层都有一定数量的神经元，不同数量的神经元组成的多层网络有很强大的非线性表达能力。以三层BPNN为例，其网络结构如下图所示。在偏差订正方面，BPNN能够捕捉到模式结果与观测之间偏差的非线性变化，得到更精准的模式订正结果。因此，我们选择了这个模型来做后续的订正。

BPNN 模型订正历史全球海温数据

starttime = time.time()
for n in range(6):
    #### 归一化
    ### 训练集
    Scaler_Input = MinMaxScaler(feature_range = (-1, 1) )
    Input_Train_Scaler[n,:,:] = Scaler_Input.fit_transform(IMFs_Historical_Compose_3D_Train[n,:,:])
    Scaler_Output = MinMaxScaler(feature_range = (-1, 1) )
    Output_Train_Scaler[n,:,:] = Scaler_Output.fit_transform(IMFs_ERSST_Compose_3D_Train[n,:,:])
    Input_Validate_Scaler[n,:,:] = Scaler_Input.transform(IMFs_Historical_Compose_3D_Validate[n,:,:])
    Input_Test_Scaler[n,:,:] = Scaler_Input.transform(IMFs_Historical_Compose_3D_Test[n,:,:])
    Input_Forecast126_Scaler[n,:,:] =  Scaler_Input.transform(IMFs_SSP126_Compose_3D[n,:,:])
    Input_Forecast245_Scaler[n,:,:] =  Scaler_Input.transform(IMFs_SSP245_Compose_3D[n,:,:])
    Input_Forecast585_Scaler[n,:,:] =  Scaler_Input.transform(IMFs_SSP585_Compose_3D[n,:,:])

    #### 训练模型，尝试不同隐层神经元
    for number in range(6,15):  # 从6到15对隐藏层神经元调试
        ### 建立模型
        Model_MLPR = MLPRegressor(solver = 'sgd',
                                  activation = 'relu',
                                  alpha = 1e-4,
                                  learning_rate_init = 0.01,
                                  max_iter = 200,
                                  hidden_layer_sizes = (number),
                                  random_state = 1  )
        ### 训练网络
        Model_MLPR.fit( Input_Train_Scaler[n,:,:], Output_Train_Scaler[n,:,:] )  # 输入: (样本条数, 输入层神经元数); 输出：(样本条数, 输出层神经元数)
        Predict_Train_Scaler = Model_MLPR.predict(Input_Train_Scaler[n, :, :])  # (训练样本年份数, 输入神经元数)
        ## 反归一化
        Predict_Train_Real = Scaler_Output.inverse_transform( Predict_Train_Scaler )   # (训练样本年份数, 输入神经元数)
        Predict_Train_Real_1D = Predict_Train_Real.reshape(-1)   # (训练样本月份数, )
        ## 计算训练RMSE
        RMSE_Train = rmse(IMFs_Historical_Compose_2D_Train[n, :], IMFs_ERSST_Compose_2D_Train[n,:] )
        RMSE_Predict_Train_All[n, number-6] = rmse(Predict_Train_Real_1D, IMFs_ERSST_Compose_2D_Train[n,:] )
        print('IMF number is: ', n + 1, ', neuron number is: ',number,  ', RMSE_Train: ', RMSE_Train, ',RMSE_Predict_Train: ', RMSE_Predict_Train_All[n, number-6])
       
        ### 验证集
        Predict_Validate_Scaler = Model_MLPR.predict(Input_Validate_Scaler[n, :, :])  # (验证样本年份数,输出层神经元数 )
        ## 反归一化
        Predict_Validate_Real = Scaler_Output.inverse_transform( Predict_Validate_Scaler )   # (验证样本年份数, 输出层神经元数)
        Predict_Validate_Real_1D = Predict_Validate_Real.reshape(-1)  # (验证样本月份数, )
        ## 计算验证集RMSE
        RMSE_Validate = rmse(IMFs_Historical_Compose_2D_Validate[n, :], IMFs_ERSST_Compose_2D_Validate[n,:] )
        RMSE_Predict_Validate_All[n, number - 6] = rmse(Predict_Validate_Real_1D, IMFs_ERSST_Compose_2D_Validate[n, :])
        print('IMF number is: ', n + 1, ', neuron number is: ',number, ', RMSE_Validate: ', RMSE_Validate, ',RMSE_Predict_Validate: ', RMSE_Predict_Validate_All[n, number-6])
        
        ### 测试集
        Predict_Test_Scaler = Model_MLPR.predict(Input_Test_Scaler[n, :, :])  # (测试样本年份数, 输出层神经元数)
        ## 反归一化
        Predict_Test_Real = Scaler_Output.inverse_transform(Predict_Test_Scaler)  # (测试样本年份数, 输出层神经元数)
        Predict_Test_Real_1D = Predict_Test_Real.reshape(-1)  # (测试样本月份数, )
        ## 计算测试集RMSE
        RMSE_Test = rmse(IMFs_Historical_Compose_2D_Test[n, :], IMFs_ERSST_Compose_2D_Test[n, :])
        RMSE_Predict_Test_All[n, number - 6] = rmse(Predict_Test_Real_1D, IMFs_ERSST_Compose_2D_Test[n, :])
        print('IMF number is: ', n + 1, ', neuron number is: ',number,  ', RMSE_Test: ', RMSE_Test, ',RMSE_Predict_Test: ', RMSE_Predict_Test_All[n, number - 6])
        
        ### SSP：未来预估订正
        Predict_SSP126_Scaler = Model_MLPR.predict(Input_Forecast126_Scaler[n, :, :])  # (未来预估年份数, 输出层神经元数)
        Predict_SSP245_Scaler = Model_MLPR.predict(Input_Forecast245_Scaler[n, :, :])  # (未来预估年份数, 输出层神经元数)
        Predict_SSP585_Scaler = Model_MLPR.predict(Input_Forecast585_Scaler[n, :, :])  # (未来预估年份数, 输出层神经元数)


    

        ## 反归一化
        Predict_SSP126_Real = Scaler_Output.inverse_transform(Predict_SSP126_Scaler)  # (未来预估年份数, 输出层神经元数)
        Predict_SSP245_Real = Scaler_Output.inverse_transform(Predict_SSP245_Scaler)  # (未来预估年份数, 输出层神经元数)
        Predict_SSP585_Real = Scaler_Output.inverse_transform(Predict_SSP585_Scaler)  # (未来预估年份数, 输出层神经元数)

        ### 找到验证集最小的RMSE及其对应的隐层神经元数，以保存所有集合的订正结果
        if number == 6:
            RMSE_Predict_Validate_Min = rmse(Predict_Validate_Real_1D, IMFs_ERSST_Compose_2D_Validate[n, :] )
            print('--------------------------------------------------------------------------------------------------------------------------------')
            print('IMF number is: ', n + 1, ', neuron number is: ',number,  ', RMSE_Predict_Validate_Min: ', RMSE_Predict_Validate_Min)
            print('--------------------------------------------------------------------------------------------------------------------------------')
            Neuron_Predict_All[n] = number
            BP_Output_Train[n,:,:] = Predict_Train_Real
            BP_Output_Validate[n,:,:] = Predict_Validate_Real
            BP_Output_Test[n,:,:] = Predict_Test_Real
            BP_Output_SSP126[n,:,:] = Predict_SSP126_Real
            BP_Output_SSP245[n,:,:] = Predict_SSP245_Real
            BP_Output_SSP585[n,:,:] = Predict_SSP585_Real
            exec('Model_MLPR_IMF_%s = Model_MLPR'%(n + 1) )
            exec("pickle.dump(Model_MLPR_IMF_%s, open(Path_SaveModel + 'Model_MLPR_IMF_%s', 'wb') )"%(n + 1, n + 1) )
        elif number >= 7:
            RMSE_Predict_Validate_Current = rmse(Predict_Validate_Real_1D, IMFs_ERSST_Compose_2D_Validate[n, :] )
            if RMSE_Predict_Validate_Current                 RMSE_Predict_Validate_Min = RMSE_Predict_Validate_Current
                print('--------------------------------------------------------------------------------------------------------------------------------')
                print('IMF number is: ', n + 1, ', neuron number is: ', number, ', RMSE_Predict_Validate_Min: ', RMSE_Predict_Validate_Min)
                print('--------------------------------------------------------------------------------------------------------------------------------')
                Neuron_Predict_All[n] = number
                BP_Output_Train[n, :, :] = Predict_Train_Real
                BP_Output_Validate[n, :, :] = Predict_Validate_Real
                BP_Output_Test[n, :, :] = Predict_Test_Real
                BP_Output_SSP126[n, :, :] = Predict_SSP126_Real
                BP_Output_SSP245[n, :, :] = Predict_SSP245_Real
                BP_Output_SSP585[n, :, :] = Predict_SSP585_Real
                exec('Model_MLPR_IMF_%s = Model_MLPR' % (n + 1))
                exec("pickle.dump(Model_MLPR_IMF_%s, open(Path_SaveModel + 'Model_MLPR_IMF_%s', 'wb') )"%(n + 1, n + 1) )


endtime = time.time()
print('spend time is ', str(endtime - starttime) + "s")

> spend time is  1.8166329860687256s

# 打印每个组合验证集修正后误差最小时对应的隐层神经元数目
print(Neuron_Predict_All)

> [10.  7. 12. 14.  7. 14.]

BPNN订正结果分析

条形统计图

未来预估的时间序列