一种气象台站测量日总辐射曝辐量预测优化方法与流程

本发明属于太阳能
技术领域：
，特别涉及了一种气象台站日总辐射曝辐量预测优化方法。
背景技术：
：太阳辐射是地球大气运动的主要能源，也是地球气候形成的最重要因子之一。多年来的实践证明，太阳辐射观测资料在对大气运动规律研究、气候预测、农作物产量评估、气候资源的开发利用等研究领域是必不可少的基础数据。日总辐射曝辐量作为太阳辐射的要素之一，其预报方法也经历了早期根据天气形势分析为主的经验方法过渡到根据历史日总辐射曝辐量记录进行的统计预测，如多元线性回归方法、逐步回归、主分量分析、奇异值分解、典型相关分析等。伴随着天气动力学及计算机技术的发展，模式预报逐渐成为过去一段时间天气预报的主要技术之一，大气模式是一套能描述气候系统中的各种物理、化学和生物过程及其相互作用的数学方程组，通过求解方程可得到温度、日总辐射曝辐量等各物理量值，但是由于天气系统的高度复杂性，模式不完备和初始场误差等因素，现有的模式预报还有改进的需要。同时，由于模式预报的复杂性，不是任何单位能够独自的进行符合本单位需要的个性化预报。随着气象信息化程度的日益提高，气象部门积累了大量的气象数据，如何充分利用这些数据，提高天气预报的精准度，提供个性化预报成为当代天气预报研究的一个重点方向。经过多年的实践与探索，以神经网络为代表的一类智能建模方法已广泛应用于天气预报中，如气温预测、降水预测等，成为模式预报的有益补充。在日总辐射曝辐量变化预测过程中，不仅预测方法关键，即使采用相同的预测方法其预测效果也不尽相同，这与具体的预测方法无关，而是天气系统的高度非线性、对初值的敏感性造成的。技术实现要素：针对上述问题，本发明提出一种气象台站日总辐射曝辐量预测优化方法，通过小波功率谱技术提取测量日总辐射曝辐量序列中的显著周期序列，分离出残差序列，由于显著周期序列规律性强，因此对显著周期序列直接建立粒子群优化的神经网络模型进行预测；而对于残差序列首先经由小波分解得到各分量序列，由于这些分量序列都归结为某个特定的频域范围，其在对应频域范围内具有近似的特征和规律，而后对各分量序列分别建立粒子群优化的神经网络模型进行预测，可以有效提高残差序列的预报精度。整体上，由于显著周期序列占原日总辐射曝辐量测量序列的比重高，规律性强，因而可以高精度的预报，而残差序列经由小波分解从而可以进一步提高对于残差序列的预报精度，从而整体提高台站测量日总辐射曝辐量的预报精度。实现上述技术目的，达到上述技术效果，本发明通过以下技术方案实现：一种气象台站测量日总辐射曝辐量预测优化方法，包括以下步骤：(1)读入经数据检验之后的气象台站测量日总辐射曝辐量时间序列，将它调整为按预报间隔要求的平均日总辐射曝辐量时间序列，然后计算出平均日总辐射曝辐量时间序列的距平序列；(2)采用多尺度小波功率谱分析方法，提取平均日总辐射曝辐量时间序列的距平序列中隐含的显著周期序列，并分离得到残差序列；(3)对平均日总辐射曝辐量时间序列的距平序列中的显著周期序列分别采用粒子群算法优化bp神经网络进行预测，获得各显著周期序列的预测结果；(4)对残差序列r经d阶小波分解得到{ad,d1,d2,…,di,...,dd}(i＝1,2,...,d)，共d+1个分量序列，对这d+1个分量序列分别采用粒子群算法优化bp神经网络进行预测，并获得相应的预测结果；(5)将平均日总辐射曝辐量时间序列的平均值与各显著周期序列、残差序列的预测结果相加，得到最终的预测结果。进一步地，所述气象台站测量日总辐射曝辐量时间序列为p＝{p(i),i＝1,2,...,n}，其中n为气象台站原始测量日总辐射曝辐量序列采样点个数；所述平均日总辐射曝辐量时间序列为p’＝{p’(j),j＝1,2,...,m}，其中m为按预报间隔要求调整后的气象台站平均测量日总辐射曝辐量序列的采样点个数；p’的平均值为令所述平均日总辐射曝辐量时间序列的距平序列为进一步地，所述显著周期序列为{p1,p2,…,pk,…,pk}，其中k为p中隐含的显著周期序列的个数，pk＝{pk(1),pk(2),…,pk(m)}，其中pk(1),pk(2),…,pk(m)分别为显著周期序列pk的值；所述残差序列为r＝p-p1-p2-…-pk。进一步地，所述采用小波功率谱分析方法，提取平均测量日总辐射曝辐量序列的距平序列中隐含的显著周期序列，具体为：利用小波功率谱方法，分析平均测量日总辐射曝辐量序列的距平序列的显著周期带，并利用小波重构方法提取各显著周期带对应的时间序列，从而获得显著周期序列。进一步地，所述步骤(3)具体为：设定预测步长为l，则各显著周期序列{p1,p2,…,pk,…,pk}的预测结果为其中进一步地，所述步骤(4)具体为：设定预测步长为l，残差序列r的小波分解分量序列{ad,d1,d2,…,di,...,dd}的预测结果为其中则残差序列r的预测结果为进一步地，所述步骤(3)中针对显著周期序列所采用的粒子群算法优化bp神经网络的具体过程为：(1)依据kolmogorov定理，建立3层bp神经网络模型，设输入层神经元个数为i，隐含层神经元个数为h，输出层神经元个数为o；其中，h＝2*i+1，o＝1；(2)确定需要优化的参数：包括：bp神经网络的输入层神经元个数i和训练集的长度l，还包括：一组对象w＝(w(1),w(2),...,w(q))，q＝i*h+h*o+h+o，其中，w(1)～w(i*h)为bp神经网络的输入层至隐含层神经元的连结权值，w(i*h+1)～w(i*h+h*o)为bp神经网络的隐含层至输出层神经元的连结权值，w(i*h+h*o+1)～w(i*h+h*o+h)为bp神经网络隐含层神经元的阈值，w(i*h+h*o+h+1)～w(i*h+h*o+h+o)为bp神经网络输出层神经元的阈值；(3)初始化种群x＝(x1,x2,...,xq1)，其中q1为粒子的总数，第i个粒子为xi＝(ii,wi,li)，粒子速度为vi＝(v_ii,v_wi,v_li)，其中ii、wi、li为参数i、w、l一组备选解；(4)对群体中的每个粒子xi＝(ii,wi,li)确定的参数，构造bp神经网络训练集的输入和输出矩阵，其中针对显著周期序列pk及bp神经网络输入层神经元个数ii首先建立矩阵z1和z2，其中：针对待优化神经网络训练集长度l，z1中最后的li列作为训练集的输入矩阵itrain，z2中最后的li列作为训练集的输出矩阵otrain；将预报步长l作为测试步长，z1中最后的l列作为测试集的输入矩阵itest，z2中最后的l列作为测试集的输出矩阵otest；根据训练集构造的bp神经网络对测试集模拟结果的误差平方和作为其适应度值，以适应度值最小为优化方向作为评价标准评判各个粒子的优劣，记录粒子xi当前个体极值为pbest(i)，取群体中pbest(i)最优的个体作为整体极值gbest；(5)群体中的每个粒子xi，分别对其位置和速度进行更新；式中：ω为惯性权重，c1、c2为加速度因子，g为当前迭代次数，而r1、r2为分布于[0,1]的随机数；(6)重新计算各个粒子此时的目标函数值，更新pbest(i)和gbest；(7)判断是否达到最大迭代次数，如满足则结束优化过程，获得经粒子群算法优化得到的参数最优值为(ibest,wbest(wbest(1),wbest(2),...,wbest(q)),lbest)，否则返回步骤(4)；(8)按ibest、wbest(wbest(1),wbest(2),...,wbest(q))、lbest构造bp神经网络训练集z3和测试集z4并初始化bp神经网络连结权值和阈值，其中：wbest(1)～wbest(i*h)为bp神经网络的输入层至隐含层神经元的连结权值的初始值，wbest(i*h+1)～wbest(i*h+h*o)为bp神经网络的隐含层至输出层神经元的连结权值的初始值，wbest(i*h+h*o+1)～wbest(i*h+h*o+h)为bp神经网络隐含层神经元的阈值的初始值，wbest(i*h+h*o+h+1)～wbest(i*h+h*o+h+o)为bp神经网络输出层神经元的阈值的初始值，就此建立起bp神经网络模型，经训练后进行迭代的l步预测，并获得对应的预测结果。进一步地，所述惯性权重ω＝0.5，加速度因子c1＝c2＝1.49445。进一步地，所述步骤(4)中针对小波分解分量序列所采用的粒子群算法优化bp神经网络的具体过程与所述步骤(3)中针对显著周期序列所采用的粒子群算法优化bp神经网络的具体过程相同。本发明的有益效果：(1)经多尺度小波功率谱分析提取的台站测量日总辐射曝辐量显著周期序列由于规律性强，因此可以高精度的进行预测，而且显著周期序列在原序列中所占比重大，因此奠定了较高精度预测的基础；剔除了显著周期序列后的残差序列一方面由于在原序列中的比重不大，另一方面采用小波分解将残差序列分解为各个频域范围内的分量序列，分量序列在各个频域范围内具有近似的特征和规律，因此可以有效提高对于残差序列的预报精度，从而整体提高台站测量日总辐射曝辐量的预报精度。(2)针对神经网络结构选择不一对于预报性能的影响，本发明提出针对提取的显著周期序列和分离的残差序列小波分解后得到的分量序列，采用粒子群算法优化方法，对于神经网络的结构参数，训练集规模进行优化，显著改善了神经网络的泛化性能，最终提高了预测精度。附图说明图1为本发明的方法流程图；图2为原始台站测量日总辐射曝辐量序列；图3为平均日总辐射曝辐量时间序列的距平序列的小波功率谱分析结果；图4为平均日总辐射曝辐量时间序列的距平序列提取的显著周期序列p1及分离的残差序列r；图5为残差序列r的3阶db3小波分解分量；图6(a)为本发明方法的一步预测结果；图6(b)为本发明方法的二步预测结果；图6(c)为本发明方法的三步预测结果；图7为平均日总辐射曝辐量时间序列的距平序列3阶db3小波分解分量；图8(a)为原始台站测量日总辐射曝辐量序列3阶db3小波分解-粒子群优化bp神经网络一步预测结果；图8(b)为原始台站测量日总辐射曝辐量序列3阶db3小波分解-粒子群优化bp神经网络二步预测结果；图8(c)为原始台站测量日总辐射曝辐量序列3阶db3小波分解-粒子群优化bp神经网络三步预测结果；图9(a)为原始台站测量日总辐射曝辐量序列arima时间序列模型一步预测结果；图9(b)为原始台站测量日总辐射曝辐量序列arima时间序列模型二步预测结果；图9(c)为原始台站测量日总辐射曝辐量序列arima时间序列模型三步预测结果。具体实施方式为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。下面结合附图对本发明的应用原理作详细的描述。如图1所示，一种气象台站测量日总辐射曝辐量预测优化方法，具体包括以下步骤：s1：读入经数据检验之后的气象台站测量日总辐射曝辐量时间序列，将它调整为按预报间隔要求的平均日总辐射曝辐量时间序列，然后计算出平均日总辐射曝辐量时间序列的距平序列；具体地：所述气象台站测量日总辐射曝辐量时间序列为p＝{p(i),i＝1,2,...,n}，其中n为气象台站原始测量日总辐射曝辐量序列采样点个数；所述平均日总辐射曝辐量时间序列为p’＝{p’(j),j＝1,2,...,m}，其中m为按预报间隔要求调整后的气象台站平均测量日总辐射曝辐量序列的采样点个数；p’的平均值为令所述平均日总辐射曝辐量时间序列的距平序列为所述显著周期序列为{p1,p2,…,pk,…,pk}，其中k为p中隐含的显著周期序列的个数，pk＝{pk(1),pk(2),…,pk(m)}，其中pk(1),pk(2),…,pk(m)分别为显著周期序列pk的值；所述残差序列为r＝p-p1-p2-…-pk。s2：采用多尺度小波功率谱分析方法，提取平均日总辐射曝辐量时间序列的距平序列中隐含的的显著周期序列，并分离得到残差序列；其中，所采用的小波功率谱分析方法的步骤为：假设一离散时间序列xn，其中n＝1,...,n，共n个采样点，采样时间间隔δt＝1，应用morlet小波变换，分析该时间序列的显著周期，提取各显著周期带对应的时间序列。2.1确定分析周期小波变换的周期tj(图3横坐标值)与小波分析中尺度参数sj有关，考虑morlet母小波中心周期特性，此处tj＝sj。尺度参数的选择为：sj＝2jδj，其中δj＝1/4，j＝0,1,...,j，共j+1个尺度，其中j的最大值不超过jmax＝4log2(n)即可，这里j＝48。2.2确定全局小波变换谱值在第n采样点，sj尺度参数对应的局地小波变换谱值wn(sj)为：其中ψ*(·)为ψ(·)的共轭函数，对第n采样点，尺度参数sj的morlet小波基小波函数为：对wn(sj)的模|wn(sj)|沿整个采样区间积分，得到尺度参数sj对应的全局小波变换谱值即：本发明采用标准化的全局小波变换谱值(图3实线)，其中σ2为xn的方差。2.3全局小波变换谱值显著性检验通常将全局小波变换谱值曲线的极大值对应的周期定为主要周期，但是否显著，要通过显著性检验。这里将以上所得全局小波变换谱值与红噪音谱值进行比较，判断其显著性，其中红噪音谱值qk表示为：其中α为xn时间序列落后一个采样点的自相关系数，k＝0,1,…,n/2。假设全局小波变换谱值为某一非周期过程谱值，则其与红噪音谱值之比遵从被自由度ν去除的分布：其中自由度γ为去相关因子，对morlet小波，γ＝2.32。在此取0.05显著性水平，当时，该全局小波变换谱值对应的周期是显著的，其中为图3虚线。2.4提取显著周期带对应的时间序列提取某特定周期带即[t1,t2]对应的时间序列x′n，由(1)知特定周期带[t1，t2]对应的尺度参数为对于morlet小波，提取尺度参数对应的时间序列是对该尺度参数带对应进行求和，即：其中，其中为wn(sj)的实部，针对morlet小波，ψ0(0)＝π-1/4，cδ＝0.776。s3：对平均日总辐射曝辐量时间序列的距平序列中的显著周期序列分别采用粒子群算法优化bp神经网络进行预测，获得各显著周期序列的预测结果；设定预测步长为l，则各显著周期序列{p1,p2,…,pk,…,pk}的预测结果为其中进一步地，所述步骤(4)具体为：设定预测步长为l，残差序列r的小波分解分量序列{ad,d1,d2,…,di,...,dd}的预测结果为其中(i＝1,2,...,d)，则残差序列r的预测结果为其中，所述步骤(3)中针对显著周期序列所采用的粒子群算法优化bp神经网络的具体过程为：3.1依据kolmogorov定理，建立3层bp神经网络模型，设输入层神经元个数为i，隐含层神经元个数为h，输出层神经元个数为o；其中，h＝2*i+1，o＝1；3.2确定需要优化的参数：包括：bp神经网络的输入层神经元个数i和训练集的长度l，还包括：一组对象w＝(w(1),w(2),...,w(q))，q＝i*h+h*o+h+o，其中，w(1)～w(i*h)为bp神经网络的输入层至隐含层神经元的连结权值，w(i*h+1)～w(i*h+h*o)为bp神经网络的隐含层至输出层神经元的连结权值，w(i*h+h*o+1)～w(i*h+h*o+h)为bp神经网络隐含层神经元的阈值，w(i*h+h*o+h+1)～w(i*h+h*o+h+o)为bp神经网络输出层神经元的阈值；3.3初始化种群x＝(x1,x2,...,xq1)，其中q1为粒子的总数，第i个粒子为xi＝(ii,wi,li)，粒子速度为vi＝(v_ii,v_wi,v_li)，其中ii、wi、li为参数i、w、l一组备选解；3.4对群体中的每个粒子xi＝(ii,wi,li)确定的参数，构造bp神经网络训练集的输入和输出矩阵，其中针对显著周期序列pk及bp神经网络输入层神经元个数ii首先建立矩阵z1和z2，其中：针对待优化神经网络训练集长度l，z1中最后的li列作为训练集的输入矩阵itrain，z2中最后的li列作为训练集的输出矩阵otrain；将预报步长l作为测试步长，z1中最后的l列作为测试集的输入矩阵itest，z2中最后的l列作为测试集的输出矩阵otest；根据训练集构造的bp神经网络对测试集模拟结果的误差平方和作为其适应度值，以适应度值最小为优化方向作为评价标准评判各个粒子的优劣，记录粒子xi当前个体极值为pbest(i)，取群体中pbest(i)最优的个体作为整体极值gbest；3.5群体中的每个粒子xi，分别对其位置和速度进行更新；式中：ω为惯性权重，c1、c2为加速度因子，g为当前迭代次数，而r1、r2为分布于[0,1]的随机数；3.6重新计算各个粒子此时的目标函数值，更新pbest(i)和gbest；3.7判断是否达到最大迭代次数，如满足则结束优化过程，获得经粒子群算法优化得到的参数最优值为(ibest,wbest(wbest(1),wbest(2),...,wbest(q)),lbest)，否则返回步骤3.4；3.8按ibest、wbest(wbest(1),wbest(2),...,wbest(q))、lbest构造bp神经网络训练集z3和测试集z4并初始化bp神经网络连结权值和阈值，其中：wbest(1)～wbest(i*h)为bp神经网络的输入层至隐含层神经元的连结权值的初始值，wbest(i*h+1)～wbest(i*h+h*o)为bp神经网络的隐含层至输出层神经元的连结权值的初始值，wbest(i*h+h*o+1)～wbest(i*h+h*o+h)为bp神经网络隐含层神经元的阈值的初始值，wbest(i*h+h*o+h+1)～wbest(i*h+h*o+h+o)为bp神经网络输出层神经元的阈值的初始值，就此建立起bp神经网络模型，经训练后进行迭代的l步预测，并获得对应的预测结果。其中，所述惯性权重ω＝0.5，加速度因子c1＝c2＝1.49445。s4：对残差序列r经d阶小波分解得到{ad,d1,d2,…,di,...,dd}(i＝1,2,...,d)，共d+1个分量序列，对这d+1个分量序列分别采用粒子群算法优化bp神经网络进行预测，并获得相应的预测结果；所述步骤s4中针对小波分解分量序列所采用的粒子群算法优化bp神经网络的具体过程与所述步骤s3中针对显著周期序列所采用的粒子群算法优化bp神经网络的具体过程相同。s5：将平均日总辐射曝辐量时间序列的平均值与各显著周期序列、残差序列的预测结果相加，得到最终的预测结果。具体测试实例：按图1所展示的流程图，取中国气象数据网下载的自1957年开始的中国气象辐射基本要素日值数据集，所有数据已经过质量检验与控制，由于本发明给出的具体测试用例也是日周期级别，因此该数据可以直接使用。本测试实例取我国编号为51628的新疆阿克苏台站的自2012年1月1日起至2016年12月31日止的日总辐射曝辐量时间序列，共2192个逐日平均日总辐射曝辐量数据，如图2所示，取2016年1月31日之前1857个数据作为已知数据集，开展对自2016年2月1日开始的为期为50个步长的1步、2步和3步预测实验，并以相对百分比误差mape为标准测试本算法的有效性：其中，y(i)和p’(i)分别为台站日总辐射曝辐量预测值和采样值，l为预测步长。平均日总辐射曝辐量时间序列的距平序列p的小波功率谱分析结果如图3所示，以5％显著性水平的红噪音检测线为阈值，该台站测量平均日总辐射曝辐量时间序列的距平序列仅具有一个以362采样点为极值点的显著周期，即年周期，取其极值点左右两侧各第一个低于红噪音检测线的周期点，组成周期带，此周期带为显著周期带，以本例为例，本例仅存在一个显著周期带即[215.3,512]，根据小波重构方法，提取此周期带对应的显著周期序列p1，并得到对应的残差序列r，由此p＝p1+r，见图4。可见，显著周期序列p1的规律性极强，可以较高精度的预测；另一方面，虽然针对残差的预测误差不可避免，但经计算，一方面残差r的能量(方差)占比p的能量(方差)为33.35％，下降明显，另一方面，残差r经3阶db3小波分解后的分量序列，如图5所示，其近似部分a3与细节部分d1～d3则呈现围绕0轴的波动，因此对于r经3阶db3小波分解后的分量序列的预测结果有望好于直接对r进行预测的结果。对p1采用基于粒子群算法优化的bp神经网络，取输入层神经元个数的范围为[5,14]，训练集的长度为[50,1800]，神经网络权值和阈值的范围为[-3,3]，粒子群种群规模是50，迭代30次。对于r则首先采用3阶db3小波分解，生成r的近似部分a3以及细节部分d1、d2和d3，之后对{a3,d1,d2,d3}分别采用基于粒子群算法优化的bp神经网络，取输入层神经元个数的范围为[5,14]，训练集的长度为[50,1800]，神经网络权值和阈值的范围为[-3,3]，粒子群种群规模是50，迭代30次。表1所示为进行3步预测时，针对显著周期序列p1和残差r的小波分解分量序列{a3,d1,d2,d3}输入层神经元个数i和训练集长度l两个参数的优化结果，由于bp神经网络权值和阈值的优化结果参数过多而不一一列出。表1本测试用例进行了总预测步长为50的1步、2步和3步预测实验，预测结果如图6(a)-(c)所示，表2为预测误差统计。表21步预测2步预测3步预测mape0.09490.12720.1281对比实验1直接对平均日总辐射曝辐量时间序列的距平序列p进行3阶db3小波分解，如图7所示，之后对小波分解分量{a3,d1,d2,d3}分别采用基于粒子群算法优化的bp神经网络，取输入层神经元个数的范围为[5,14]，训练集的长度为[50,1800]，神经网络权值和阈值的范围为[-3,3]，粒子群种群规模是50，迭代30次。表3所示为进行3步预测时，针对小波分解分量{a3,d1,d2,d3}输入层神经元个数i和训练集长度l两个参数的优化结果，由于bp神经网络权值和阈值的优化结果参数过多而不一一列出。表3本测试用例进行了总预测步长为50的1步、2步和3步预测实验，预测结果如图8(a)-(c)所示，表4为预测误差统计，可见其1～3步的mape平均误差较表2提高了33.15％，说明首先提取显著周期序列，之后对残差序列进行小波分解的策略，可以有效的提高预测精度。表41步预测2步预测3步预测mape0.10790.15720.2012对比实验2台站原始测量日总辐射曝辐量序列建立差分自回归移动平均模型(autoregressiveintegratedmovingaveragemodel，arima)模型。选取预测点前100个采样数据点，通过aic准则定阶法确定arima模型的结构，同样的，对比实验2进行了总预测步长为50的1步、2步和3步预测实验，预测结果如图9(a)-(c)所示，表5为预测误差统计，可见其1～3步的mape平均误差较表2提高了62.91％。表5对原始序列p’建立的arima时间序列模型预测误差1步预测2步预测3步预测mape0.12190.19820.2504以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。当前第1页12