CN115146526B - 一种非平稳过程粉尘浓度预测方法 - Google Patents

一种非平稳过程粉尘浓度预测方法 Download PDF

Info

Publication number
CN115146526B
CN115146526B CN202210430150.0A CN202210430150A CN115146526B CN 115146526 B CN115146526 B CN 115146526B CN 202210430150 A CN202210430150 A CN 202210430150A CN 115146526 B CN115146526 B CN 115146526B
Authority
CN
China
Prior art keywords
trend
time
dust concentration
predict
resid
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210430150.0A
Other languages
English (en)
Other versions
CN115146526A (zh
Inventor
邓勤
吴付祥
惠立锋
赵政
张强
陈建阁
晏丹
李征真
杨亚会
李彦筑
王宇廷
罗小博
焦敏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
CCTEG Chongqing Research Institute Co Ltd
Original Assignee
CCTEG Chongqing Research Institute Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by CCTEG Chongqing Research Institute Co Ltd filed Critical CCTEG Chongqing Research Institute Co Ltd
Priority to CN202210430150.0A priority Critical patent/CN115146526B/zh
Publication of CN115146526A publication Critical patent/CN115146526A/zh
Application granted granted Critical
Publication of CN115146526B publication Critical patent/CN115146526B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
    • G01N15/06—Investigating concentration of particle suspensions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10—Noise analysis or noise optimisation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Evolutionary Computation (AREA)
  • Databases & Information Systems (AREA)
  • Chemical & Material Sciences (AREA)
  • Geometry (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Medical Informatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Dispersion Chemistry (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Operations Research (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及一种非平稳过程粉尘浓度预测方法，属于矿山粉尘预测领域。该方法包括以下步骤：S1：记录一段时间内粉尘浓度随时间的变化值C(t)；S2：利用变分模式分解法对C(t)进行分解，得到粉尘浓度随时间沉积趋势值、粉尘浓度随矿井工作状态周期变化值及具有平稳过程的粉尘浓度随时间波动值；S3：分别对粉尘浓度随时间沉积趋势值、粉尘浓度随矿井工作状态周期变化值及具有平稳过程的粉尘浓度随时间波动值进行预测，最终预测粉尘浓度。

Description

一种非平稳过程粉尘浓度预测方法
技术领域
本发明属于矿山粉尘预测领域，涉及一种非平稳过程粉尘浓度预测方法。
背景技术
粉尘主动控制技术是矿井职业危害防治的重要组成部分，而科学有依据的粉尘主动控制技术需依赖于具有高可靠性的粉尘浓度预测技术。
矿井粉尘浓度受矿井工作状态、生产工艺参数、粉尘颗粒沉降特性、接尘点位置、温度、湿度等多方面因素影响，考虑到多参数预测模型计算量庞大且过于复杂难以建立，目前主流的预测方法为仅将时间变量作为单变量来对粉尘浓度进行预测，然而常用的基于概率统计或时间预测序列的预测方法均要求用于预测的历史数据为统计学意义上的平稳过程，对于矿井粉尘而言，由于受矿井生产机械工作状态及粉尘沉降特性的影响，其浓度值在长时间范围内具有趋势性及部分周期性，呈现为统计学意义上的非平稳过程，因此常用的时序预测方法难以满足矿井粉尘预测需求。
发明内容
有鉴于此，本发明的目的在于提供一种非平稳过程粉尘浓度预测方法，来解决传统预测方法无法准确预测矿井粉尘浓度的问题；具有以下优点：模型对历史数据拟合精度高，后检验差比值小于0.35，模型预测精度高，相对误差小于5％。
为达到上述目的，本发明提供如下技术方案：
一种非平稳过程粉尘浓度预测方法，该方法包括以下步骤：
S1：记录一段时间内粉尘浓度随时间的变化值C(t),t＝t1,t2...tn；
S2：利用变分模式分解法对C(t)进行分解，得到粉尘浓度随时间沉积趋势值Ctrend(t)、粉尘浓度随矿井工作状态周期变化值Cseasonal(t)及具有平稳过程的粉尘浓度随时间波动值Cresid(t)，粉尘浓度C(t)＝Ctrend(t)+Cseasonal(t)+Cresid(t)；
S3：分别对Ctrend(t)、Cseasonal(t)、Cresid(t)进行预测，最终预测粉尘浓度为Cpredict(t')＝Ctrend-predict(t')+Cseasonal-predict(t')+Cresid-predict(t'),t'＝tn+1,tn+2,...,tn+m。
可选的，所述S3具体为：
利用灰色预测法对粉尘浓度随时间沉积趋势值Ctrend(t)进行预测，得到Ctrend-predict(t')；
(1)建立灰色模型；
Ctrend-predict(t)的累加序列为
Figure GDA0004209626830000021
根据灰色模型理论，Ctrend-predict (1)(t)满足方程微分方程：
Figure GDA0004209626830000022
其中，θ为发展灰数，χ为内生控制灰数；
对微分方程求解得：
Figure GDA0004209626830000023
对Ctrend-predict (1)(t)进行逆生还原得：
Ctrend-predict(tn+m)＝Ctrend (1)(tn+m)-Ctrend (1)(tn+m-1)
(2)利用已知数据求解θ，χ；
由于t为离散时间序列，故
Figure GDA0004209626830000024
写为：
Figure GDA0004209626830000025
为提高模型拟合与预测精度，θ，χ的取值原则为：令Ctrend-predict(t)的值在已知时间序列中逼近Ctrend(t)；在求解过程中令
Figure GDA0004209626830000026
同时将项Ctrend-predict (1)(t)修正为均值生成序列：
Figure GDA0004209626830000027
则微分方程写为：
Ctrend(t)+θZtrend (1)(t)＝χ
利用最小二乘法求得：
Figure GDA0004209626830000028
其中：
Figure GDA0004209626830000031
可选的，所述S3具体为：
利用三次指数平滑算法，即Holter-Winter法对粉尘浓度随矿井工作状态周期变化值Cseasonal(t)进行预测，得到Cseasonal-predict(t')；
Cseasonal-predict(tn+m)＝(1-κ)[A(tn)+B(tn)(tn+m-tn)+S(tn+m-T+1+mod((tn+m-tn)/T))]其中T为Cseasonal(t)变化周期，A(tn)为一次平滑值，B(tn)为二次平滑值，S(tn)为三次平滑值，计算公式如下：
A(tn)＝α'(Cseasonal(tn)-S(tn-T))+(1-α')(A(tn-1)+B(tn-1))
B(tn)＝β(A(tn)-A(tn-1))+(1-β)B(tn-1)
S(tn)＝γ(Cseasonal(tn)-A(tn-1)-B(tn-1))+(1-γ)S(tn-T)
其中α'、β、γ为平滑系数，利用K-Fold交叉验证法以RMSE最小化准则求得；
(1-κ)为引入的开关变量，用于衡量矿井工作状态，κ∈[0,1]，其取值与矿井工作强度呈反比，若预测时间段矿井工作状态与已知时间段相同，则κ＝0，若预测时间段矿井处于停工状态，则κ＝1。
可选的，所述S3具体为：
利用ARMA法对对平稳过程序列Cresid(t)项进行建模预测，得到Cresid-predict(t')；
利用变分模式分解VMD得到Ctrend(t)、Cseasonal(t)、Cresid(t)；
VMD利用迭代搜索和求解变分模型最优解来将原始复杂输入信号分解为K个中心频率为ωk的调幅调频分量信号uk(t)，称为本征模态函数IMF，每个模态在解调成基带后是平滑的，适用于时序预测方法；
令K＝3，得到IMF1、IMF2、IMF3，其分别对应Ctrend(t)、Cseasonal(t)、Cresid(t)；
各个IMF求解过程如下：
(1)构造变分问题
为保证分解序列为具有中心频率的有限带宽模态分量，同时各个模态估计带宽之和最小，约束条件为所有模态之和与原始信号相等，则VMD约束变分模型如下：
Figure GDA0004209626830000041
Figure GDA0004209626830000042
(2)求解变分问题，得到IMF
为求解此约束最优化问题，利用二次惩罚项和拉格朗日乘子法优势，通过引入增广Lagrangian函数将其转变为非约束变分问题求解：
Figure GDA0004209626830000043
其中α为二次惩罚因子，λ为拉格朗日乘子；
利用交替方向乘子法ADMM求解该非约束变分问题：
①设迭代次数为n’，初始化{uk 1},{ωk 1},λ1,n'＝0，通过傅里叶变换计算得到各个参量的频域值，计为
Figure GDA0004209626830000044
②n'＝n'+1,对k＝1:3，更新
Figure GDA0004209626830000045
ωk；
Figure GDA0004209626830000046
Figure GDA0004209626830000047
③更新
Figure GDA0004209626830000048
Figure GDA0004209626830000049
/>
其中τ表示噪声容量，当信号含有强噪声时，设定τ＝0达到更好的去噪效果；
④计算是否满足迭代约束条件：
Figure GDA0004209626830000051
其中ε为设定精度；
若满足，则停止迭代，3个IMF即为输出结果，即IMF1＝Ctrend(t),IMF2＝Cseasonal(t),IMF3＝Cresid(t)，否则返回步骤②；
利用自回归滑动平均ARMA模型对Cresid(t)项进行预测；
ARMA模型为自回归模型AR与移动平均模型MA的结合，通过记录某个或某组变量在一系列时刻的取值，得到离散数字组成的序列集合，称之为时间序列，预测过程如下：
(1)记录被观测系统历史数据，Cresid(t)；
(2)建立预测模型
Figure GDA0004209626830000052
其中μ为常数，p为模型采用时序数据本身的滞后数，称为AR项，ai为Cresid(t)序列i阶AR系数，Cresid(t-i)为Cresid(t)滞后i阶序列，ε(t)为均值为零的独立同分布白噪声，q为模型采用预测误差的滞后阶数，也称为MA项，bi为序列ε(t)i阶MA系数，ε(t-i)为ε(t)滞后i阶序列；
(3)模型定阶，求解p、q；
①计算时序数据的偏自相关系数θi，绘制其以滞后阶数i为变量的曲线图，观测曲线与置信区间
Figure GDA0004209626830000053
的关系，确定p阶数上限P＝pmax；
②计算时序数据的自相关系数ρi，绘制其以滞后阶数i为变量的曲线图，观测曲线与置信区间
Figure GDA0004209626830000054
确定q阶数上限Q＝qmax；
③对每个(p,q),0≤p≤P,0≤q≤Q，计算其对应的AIC及BIC函数值，取令AIC及BIC最小的(p,q)，定位模型阶数：
Figure GDA0004209626830000055
Figure GDA0004209626830000056
(4)利用模型进行数据预测
将待预测的时间序列t'＝tn+1,tn+2,...,tn+m代入ARMA模型中，的到Cresid-predict(t')。
本发明的有益效果在于：根据矿井粉尘特点将矿井粉尘浓度分解为三个分量，同时通过分析三个分量变化特征选择合适的浓度预测算法分别对其浓度值进行预测，本方法运算速度快，模型拟合效果极佳(后检验差比值<0.35)，且预测精度高(最大相对误差<5％)，可实现矿井粉尘的高精度预测。
本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。
附图说明
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：
图1为本发明流程图；
图2为粉尘浓度随时间变化曲线图；
图3为VMD对粉尘浓度随时间变化曲线进行分解的对应图形；
图4为原始值与利用灰色预测法对Ctrend进行预测的预测值比较曲线；
图5为原始值与利用Holter-Winters算法对Cseasonal进行预测的预测值比较曲线；
图6为原始值与利用ARMA算法对Cresid进行预测的预测值比较曲线；
图7为原始值与利用粉尘浓度已知模型拟合预测的预测值比较曲线。
具体实施方式
以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。
其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本发明的限制；为了更好地说明本发明的实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。
本发明实施例的附图中相同或相似的标号对应相同或相似的部件；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”、“前”、“后”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本发明的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。
请参阅图1，为一种非平稳过程粉尘浓度预测方法：
S1：记录一段时间内粉尘浓度随时间的变化值C(t),t＝t1,t2...tn；
S2：利用变分模式分解法对C(t)进行分解，得到粉尘浓度随时间沉积趋势值Ctrend(t)、粉尘浓度随矿井工作状态周期变化值Cseasonal(t)及具有平稳过程的粉尘浓度随时间波动值Cresid(t)，粉尘浓度C(t)＝Ctrend(t)+Cseasonal(t)+Cresid(t)；
S3：分别对Ctrend(t)、Cseasonal(t)、Cresid(t)进行预测，最终预测粉尘浓度为Cpredict(t')＝Ctrend-predict(t')+Cseasonal-predict(t')+Cresid-predict(t'),t'＝tn+1,tn+2,...,tn+m；
S31：利用灰色预测法对粉尘浓度随时间沉积趋势值Ctrend(t)进行预测，得到Ctrend-predict(t')；
(1)建立灰色模型
Ctrend-predict(t)的累加序列为
Figure GDA0004209626830000071
根据灰色模型理论，Ctrend-predict (1)(t)满足方程微分方程：
Figure GDA0004209626830000072
其中θ为发展灰数，χ为内生控制灰数
对微分方程求解可得：
Figure GDA0004209626830000073
对Ctrend-predict (1)(t)进行逆生还原可得：
Ctrend-predict(tn+m)＝Ctrend (1)(tn+m)-Ctrend (1)(tn+m-1)
(2)利用已知数据求解θ,χ
由于t为离散时间序列，故
Figure GDA0004209626830000081
可写为：
Figure GDA0004209626830000082
为提高模型拟合与预测精度，θ,χ的取值原则为：令Ctrend-predict(t)的值在已知时间序列中尽可能逼近Ctrend(t)。故在求解过程中令
Figure GDA0004209626830000083
同时将项Ctrend-predict (1)(t)修正为均值生成序列：
Figure GDA0004209626830000084
则微分方程可写为：
Ctrend(t)+θZtrend (1)(t)＝χ
利用最小二乘法可求得：
Figure GDA0004209626830000085
其中：
Figure GDA0004209626830000086
S32：利用三次指数平滑算法(Holter-Winter法)对粉尘浓度随矿井工作状态周期变化值Cseasonal(t)进行预测，得到Cseasonal-predict(t')；
Cseasonal-predict(tn+m)＝(1-κ)[A(tn)+B(tn)(tn+m-tn)+S(tn+m-T+1+mod((tn+m-tn)/T))]其中T为Cseasonal(t)变化周期，A(tn)为一次平滑值，B(tn)为二次平滑值，S(tn)为三次平滑值，计算公式如下：
A(tn)＝α'(Cseasonal(tn)-S(tn-T))+(1-α')(A(tn-1)+B(tn-1))
B(tn)＝β(A(tn)-A(tn-1))+(1-β)B(tn-1)
S(tn)＝γ(Cseasonal(tn)-A(tn-1)-B(tn-1))+(1-γ)S(tn-T)
其中α'、β、γ为平滑系数，利用K-Fold交叉验证法以RMSE最小化准则求得。
(1-κ)为本发明引入的开关变量，用于衡量矿井工作状态，κ∈[0,1]，其取值与矿井工作强度呈反比，若预测时间段矿井工作状态与已知时间段相同，则κ＝0，若预测时间段矿井处于停工状态，则κ＝1。
S33：利用ARMA法对对平稳过程序列Cresid(t)项进行建模预测，得到Cresid-predict(t')；
利用变分模式分解(VMD)得到Ctrend(t)、Cseasonal(t)、Cresid(t)；
VMD是由Konstantin Dragomiretskiy于2014年提出的一种多尺度频率分解方法，属于完全非递归模型，利用迭代搜索和求解变分模型最优解来将原始复杂输入信号分解为K个中心频率为ωk的调幅调频分量信号uk(t)，并将之称为本征模态函数(Intrinsic ModeFunction，IMF)，每个模态在解调成基带后是平滑的，适用于时序预测方法。
本发明中令K＝3，得到IMF1、IMF2、IMF3，其分别对应Ctrend(t)、Cseasonal(t)、Cresid(t)；
各个IMF求解过程如下：
(1)构造变分问题
为保证分解序列为具有中心频率的有限带宽模态分量，同时各个模态估计带宽之和最小，约束条件为所有模态之和与原始信号相等，则VMD约束变分模型如下：
Figure GDA0004209626830000091
Figure GDA0004209626830000092
(2)求解变分问题，得到IMF
为求解此约束最优化问题，可利用二次惩罚项和拉格朗日乘子法优势，通过引入增广Lagrangian函数将其转变为非约束变分问题求解
Figure GDA0004209626830000093
其中α为二次惩罚因子，λ为拉格朗日乘子
利用交替方向乘子法(Alternating Direction Multiplier Method,ADMM)求解该非约束变分问题：
①设迭代次数为n’，初始化{uk 1},{ωk 1},λ1,n'＝0，通过傅里叶变换计算得到各个参量的频域值，计为
Figure GDA0004209626830000101
②n'＝n'+1,对k＝1:3，更新
Figure GDA0004209626830000102
ωk；
Figure GDA0004209626830000103
Figure GDA0004209626830000104
③更新
Figure GDA0004209626830000105
Figure GDA0004209626830000106
其中τ表示噪声容量，当信号含有强噪声时，可设定τ＝0达到更好的去噪效果；
④计算是否满足迭代约束条件：
Figure GDA0004209626830000107
其中ε为设定精度；
若满足，则停止迭代，此时的3个IMF即为输出结果，即IMF1＝Ctrend(t),IMF2＝Cseasonal(t),IMF3＝Cresid(t)，否则返回步骤②。
S34：利用自回归滑动平均(ARMA)模型对Cresid(t)项进行预测
ARMA模型为自回归模型(AR)与移动平均模型(MA)的结合，此模型通过记录某个或某组变量在一系列时刻的取值，得到离散数字组成的序列集合，称之为时间序列，对于具有统计学意义上平稳过程的时间序列，该方法认为可通过寻找历史数据间的相关性来预测未来数据，预测过程如下：
(1)记录被观测系统历史数据，Cresid(t)；
(2)建立预测模型
Figure GDA0004209626830000111
其中μ为常数，p为模型采用时序数据本身的滞后数，也称为AR项，ai为Cresid(t)序列i阶AR系数，Cresid(t-i)为Cresid(t)滞后i阶序列，ε(t)为均值为零的独立同分布白噪声，q为模型采用预测误差的滞后阶数，也称为MA项，bi为序列ε(t)i阶MA系数，ε(t-i)为ε(t)滞后i阶序列
(3)模型定阶(求解p、q)
①计算时序数据的偏自相关系数θi，绘制其以滞后阶数i为变量的曲线图，观测曲线与置信区间
Figure GDA0004209626830000112
关系，确定p阶数上限P＝pmax；
②计算时序数据的自相关系数ρi，绘制其以滞后阶数i为变量的曲线图，观测曲线与置信区间
Figure GDA0004209626830000113
确定q阶数上限Q＝qmax；
③对每个(p,q),0≤p≤P,0≤q≤Q，计算其对应的AIC及BIC函数值，取令AIC及BIC最小的(p,q)，定位模型阶数：
Figure GDA0004209626830000114
Figure GDA0004209626830000115
(4)利用模型进行数据预测
将待预测的时间序列t'＝tn+1,tn+2,...,tn+m代入ARMA模型中，即可得到Cresid-predict(t')。
实施例
某个时间段粉尘浓度值如表1所示。
表1某时间段粉尘浓度值
Figure GDA0004209626830000116
/>
Figure GDA0004209626830000121
粉尘浓度随时间变化曲线如图2所示。
利用VMD对粉尘浓度随时间变化曲线进行分解，分别得到Ctrend、Cseasonal、Cresid分量，数值如表2所示。
表2VMD分解的分量数值
Figure GDA0004209626830000131
Figure GDA0004209626830000141
对应图形如图3所示。
令已知时序数据的前40个数据为模型训练数据集，后8个数据为模型预测验证数据集：(1)利用灰色预测法对Ctrend进行预测
①按照专利所述方法，代入训练数据集求解灰色算法模型
②利用模型对训练集及验证集进行计算，可得如表3所示的训练集。
表3利用灰色预测法对Ctrend进行预测的训练集
Figure GDA0004209626830000142
Figure GDA0004209626830000151
预测验证的数据集如表4所示。
表4预测验证数据集
时间 灰色算法预测值
t41 3.18998842500000
t42 3.18968076400000
t43 3.18933234700000
t44 3.18899278300000
t45 3.18880937200000
t46 3.18873627700000
t47 3.18876964100000
t48 3.18870610600000
原始值与预测值比较曲线如图4所示。
(2)利用Holter-Winters算法对Cseasonal进行预测
①按照专利所述方法，代入训练数据集求解Holter-Winters指数平滑模型(设κ＝0)
②利用模型对训练集及验证集进行计算，得如表5所示的训练集。
表5利用Holter-Winters算法对Cseasonal进行预测的训练集
Figure GDA0004209626830000152
/>
Figure GDA0004209626830000161
预测验证的数据集如表6所示。
表6预测验证数据集
时间 Holter-Winters算法预测值
t41 0.0101691043201441
t42 0.00668153143124810
t43 -0.00416581862971496
t44 -0.0111500733509331
t45 -0.00640631937051964
t46 0.00511796829518263
t47 0.0111122985877850
t48 0.000497705096564164
t50 0.00359867890527794
原始值与预测值比较曲线如图5所示。
(3)利用ARMA算法对Cresid进行预测
①按照专利所述方法，代入训练数据集求解ARMA模型
②利用模型对训练集及验证集进行计算，得如表7所示的训练集。
表7利用ARMA算法对Cresid进行预测的训练集
时间 ARMA算法对已知数据模拟值
t1 0.0217473820000000
t2 -8.57000000000000e-06
t3 0.0177414250000000
t4 0.00825242600000000
t5 0.00354036800000000
t6 0.0187737420000000
t7 0.0305245290000000
t8 -0.0598740950000000
t9 0.0138302830000000
t10 0.0280361370000000
t11 -0.0212464490000000
t12 -0.0209286770000000
t13 0.0596479870000000
t14 -0.0525651530000000
t15 0.00406057400000000
t16 0.0257200390000000
t17 -0.0233578490000000
t18 -0.0244302080000000
t19 0.0438784570000000
t20 -0.0257668780000000
t21 -0.0179782220000000
t22 0.0186909910000000
t23 -0.0286568140000000
t24 0.00266062000000000
t25 0.0287541450000000
t26 -0.0180171090000000
t27 -0.0115223540000000
t28 0.0301861380000000
t29 -0.0269777950000000
t30 0.0101909780000000
t31 0.0132691450000000
t32 -0.0132522260000000
t33 0.00941927700000000
t34 0.0207833810000000
t35 -0.00278769400000000
t36 -0.0147939670000000
t37 0.0301757880000000
t38 -0.0218481970000000
t39 0.0107016660000000
t40 -0.0119526980000000
预测验证的数据集如表8所示。
表8预测验证数据集
Figure GDA0004209626830000171
Figure GDA0004209626830000181
原始值与预测值比较曲线如图6所示。
(4)C_predict＝Ctrend-predict+Cseasonal-predict+Cresid-predict
故粉尘浓度已知模型拟合数据与预测数据如表9所示。
表9粉尘浓度已知模型拟合预测的训练集
Figure GDA0004209626830000182
/>
Figure GDA0004209626830000191
预测验证的数据集如表10所示。
表10预测验证数据集
时间 Holter-Winters算法预测值
t41 3.18503213029437
t42 3.19903406805183
t43 3.19776935969687
t44 3.15647385773402
t45 3.18231660752882
t46 3.19233195905763
t47 3.20022339724752
t48 3.18653277570960
t50 3.18503213029437
原始值与预测值比较曲线如图7所示。
经计算，本发明模型对训练数据集拟合均方误差MSE＝7.0466e-05，后检验比值
Figure GDA0004209626830000192
(W<0.35，一级，拟合效果非常好)，预测数据均方误差9.7917e-06，最大相对误差为0.79％，预测效果极佳。
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (3)

1.一种非平稳过程粉尘浓度预测方法，其特征在于：该方法包括以下步骤：
S1：记录一段时间内粉尘浓度随时间变化值C(t),t＝t1,t2...tn；n为记录阶段包含时刻总数量，tn为第n个时刻对应时间；
S2：利用变分模式分解法对C(t)进行分解，得到粉尘浓度随时间沉积趋势值Ctrend(t)、粉尘浓度随矿井工作状态周期变化值Cseasonal(t)及具有平稳过程的粉尘浓度随时间波动值Cresid(t)，粉尘浓度随时间变化值C(t)＝Ctrend(t)+Cseasonal(t)+Cresid(t)；
S3：分别对Ctrend(t)、Cseasonal(t)、Cresid(t)进行预测，最终预测粉尘浓度为Cpredict(t')＝Ctrend-predict(t')+Cseasonal-predict(t')+Cresid-predict(t'),t'＝tn+1,tn+2,...,tn+m；t’为需进行粉尘浓度预测的时刻，m为需要向后预测的时间段内包含时刻总数量，tn+m为第n+m个时刻对应时间；
所述S3具体为：
利用三次指数平滑算法，即Holter-Winter法对粉尘浓度随矿井工作状态周期变化值Cseasonal(t)进行预测，得到Cseasonal-predict(t')；
Cseasonal-predict(tn+m)＝(1-κ)[A(tn)+B(tn)(tn+m-tn)+S(tn+m-T+1+mod((tn+m-tn)/T))]其中T为Cseasonal(t)变化周期，A(tn)为一次平滑值，B(tn)为二次平滑值，S(tn)为三次平滑值，计算公式如下：
A(tn)＝α'(Cseasonal(tn)-S(tn-T))+(1-α')(A(tn-1)+B(tn-1))
B(tn)＝β(A(tn)-A(tn-1))+(1-β)B(tn-1)
S(tn)＝γ(Cseasonal(tn)-A(tn-1)-B(tn-1))+(1-γ)S(tn-T)
其中α'、β、γ为平滑系数，利用K-Fold交叉验证法以RMSE最小化准则求得；
(1-κ)为引入的开关变量，用于衡量矿井工作状态，κ∈[0,1]，其取值与矿井工作强度呈反比，若预测时间段矿井工作状态与已知时间段相同，则κ＝0，若预测时间段矿井处于停工状态，则κ＝1。
2.根据权利要求1所述的一种非平稳过程粉尘浓度预测方法，其特征在于：所述S3具体为：
利用灰色预测法对粉尘浓度随时间沉积趋势值Ctrend(t)进行预测，得到Ctrend-predict(t')；
(1)建立灰色模型；
Ctrend-predict(t)的累加序列为
Figure FDA0004209626820000021
tk表示第k个时刻；
根据灰色模型理论，Ctrend-predict (1)(t)满足方程微分方程：
Figure FDA0004209626820000022
其中，θ为发展灰数，χ为内生控制灰数；
对微分方程求解得：
Figure FDA0004209626820000023
对Ctrend-predict (1)(t)进行逆生还原得：
Ctrend-predict(tn+m)＝Ctrend (1)(tn+m)-Ctrend (1)(tn+m-1)
(2)利用已知数据求解θ和χ；
由于t为离散时间序列，故
Figure FDA0004209626820000024
写为：
Figure FDA0004209626820000025
为提高模型拟合与预测精度，θ,χ的取值原则为：令Ctrend-predict(t)的值在已知时间序列中逼近Ctrend(t)；在求解过程中令
Figure FDA0004209626820000026
同时将项Ctrend-predict (1)(t)修正为均值生成序列：
Figure FDA0004209626820000027
则微分方程写为：
Ctrend(t)+θZtrend (1)(t)＝χ
利用最小二乘法求得：
Figure FDA0004209626820000028
其中：
Figure FDA0004209626820000031
3.根据权利要求1所述的一种非平稳过程粉尘浓度预测方法，其特征在于：所述S3具体为：
利用ARMA法对对平稳过程序列Cresid(t)项进行建模预测，得到Cresid-predict(t')；
利用变分模式分解VMD得到Ctrend(t)、Cseasonal(t)、Cresid(t)；
VMD利用迭代搜索和求解变分模型最优解来将原始复杂输入信号分解为K个中心频率为ωk的调幅调频分量信号uk(t)，称为本征模态函数IMF，每个模态在解调成基带后是平滑的，适用于时序预测方法；
令K＝3，得到IMF1、IMF2、IMF3，其分别对应Ctrend(t)、Cseasonal(t)、Cresid(t)；
各个IMF求解过程如下：
(1)构造变分问题
为保证分解序列为具有中心频率的有限带宽模态分量，同时各个模态估计带宽之和最小，约束条件为所有模态之和与原始信号相等，则VMD约束变分模型如下：
Figure FDA0004209626820000032
/>
Figure FDA0004209626820000033
其中uk为第k个IMF序列值，ωk为第k个IMF对应的中心频率，j为虚部符号；
(2)求解变分问题，得到IMF
为求解此约束最优化问题，利用二次惩罚项和拉格朗日乘子法优势，通过引入增广Lagrangian函数将其转变为非约束变分问题求解：
Figure FDA0004209626820000034
其中L为Lagrangian函数符号，α为二次惩罚因子，λ为拉格朗日乘子；
利用交替方向乘子法ADMM求解该非约束变分问题：
①设迭代次数为n’，初始化{uk 1},{ωk 1},λ1,n'＝0，通过傅里叶变换计算得到各个参量的频域值，计为
Figure FDA0004209626820000041
②n'＝n'+1,对k＝1:3，更新
Figure FDA0004209626820000042
ωk；
Figure FDA0004209626820000043
Figure FDA0004209626820000044
③更新
Figure FDA0004209626820000045
Figure FDA0004209626820000046
其中τ表示噪声容量，当信号含有强噪声时，设定τ＝0达到更好的去噪效果；
④计算是否满足迭代约束条件：
Figure FDA0004209626820000047
其中ε为设定精度；
若满足，则停止迭代，3个IMF即为输出结果，即IMF1＝Ctrend(t),IMF2＝Cseasonal(t),IMF3＝Cresid(t)，否则返回步骤②；
利用自回归滑动平均ARMA模型对Cresid(t)项进行预测；
ARMA模型为自回归模型AR与移动平均模型MA的结合，通过记录某个或某组变量在一系列时刻的取值，得到离散数字组成的序列集合，称之为时间序列，预测过程如下：
(1)记录被观测系统历史数据，Cresid(t)；
(2)建立预测模型
Figure FDA0004209626820000048
其中μ为常数，p为模型采用时序数据本身的滞后数，称为AR项，ai为Cresid(t)序列i阶AR系数，Cresid(t-i)为Cresid(t)滞后i阶序列，ε(t)为均值为零的独立同分布白噪声，q为模型采用预测误差的滞后阶数，也称为MA项，bi为序列ε(t)i阶MA系数，ε(t-i)为ε(t)滞后i阶序列；
(3)模型定阶，求解p、q；
①计算时序数据的偏自相关系数θi，绘制其以滞后阶数i为变量的曲线图，观测曲线与置信区间
Figure FDA0004209626820000051
的关系，确定p阶数上限P＝pmax；
②计算时序数据的自相关系数ρi，绘制其以滞后阶数i为变量的曲线图，观测曲线与置信区间
Figure FDA0004209626820000052
确定q阶数上限Q＝qmax；
⑤对每个(p,q),0≤p≤P,0≤q≤Q，计算其对应的赤池信息量AIC及贝叶斯信息量BIC，取令AIC及BIC最小的(p,q)，定位模型阶数：
Figure FDA0004209626820000053
Figure FDA0004209626820000054
(4)利用模型进行数据预测
将待预测的时间序列t'＝tn+1,tn+2,...,tn+m代入ARMA模型中，的到Cresid-predict(t')。
CN202210430150.0A 2022-04-22 2022-04-22 一种非平稳过程粉尘浓度预测方法 Active CN115146526B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210430150.0A CN115146526B (zh) 2022-04-22 2022-04-22 一种非平稳过程粉尘浓度预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210430150.0A CN115146526B (zh) 2022-04-22 2022-04-22 一种非平稳过程粉尘浓度预测方法

Publications (2)

Publication Number Publication Date
CN115146526A CN115146526A (zh) 2022-10-04
CN115146526B true CN115146526B (zh) 2023-06-06

Family

ID=83407224

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210430150.0A Active CN115146526B (zh) 2022-04-22 2022-04-22 一种非平稳过程粉尘浓度预测方法

Country Status (1)

Country Link
CN (1) CN115146526B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117196109B (zh) * 2023-09-15 2024-04-05 中煤科工集团重庆研究院有限公司 一种基于多源信息融合的煤矿井下粉尘浓度预测修正方法
CN117732162A (zh) * 2024-01-24 2024-03-22 江苏海洋大学 一种煤矿矿井粉尘的过滤装置及其过滤方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109063284A (zh) * 2018-07-17 2018-12-21 西安科技大学 煤矿综掘面风筒出风口参数变化下的粉尘浓度预测方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110793896B (zh) * 2019-12-03 2022-04-08 承德石油高等专科学校 一种尾气中粉尘浓度短期预测方法
CN113834760B (zh) * 2021-10-21 2022-04-26 中国矿业大学 实时个体工作粉尘浓度暴露空间规律监测预警系统及方法
CN114166773B (zh) * 2021-12-08 2023-07-25 中煤科工集团重庆研究院有限公司 一种基于粒子群优化-支持向量机的NOx测量方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109063284A (zh) * 2018-07-17 2018-12-21 西安科技大学 煤矿综掘面风筒出风口参数变化下的粉尘浓度预测方法

Also Published As

Publication number Publication date
CN115146526A (zh) 2022-10-04

Similar Documents

Publication Publication Date Title
CN115146526B (zh) 一种非平稳过程粉尘浓度预测方法
CN110334875B (zh) 计及评估指标冲突的风电功率组合概率预测方法
CN112036084B (zh) 一种相似产品寿命迁移筛选方法和系统
CN110309603B (zh) 一种基于风速特性的短期风速预测方法及系统
US20050182624A1 (en) Method and apparatus for constructing a speech filter using estimates of clean speech and noise
CN110413754B (zh) 对话(中)奖励评估和对话方法、介质、装置和计算设备
CN106709588B (zh) 预测模型构建方法和设备以及实时预测方法和设备
Liu Instantaneous frequency tracking under model uncertainty via dynamic model averaging and particle filtering
CN109827579B (zh) 一种组合定位中滤波模型实时校正的方法和系统
CN111553513A (zh) 一种基于二次分解与回声状态网络的中长期径流预测方法
CN115601182A (zh) 一种基于改进型XGBoost类方法的数据分析方法、定价方法以及相关设备
Jia et al. Federated domain adaptation for asr with full self-supervision
CN109065176B (zh) 一种血糖预测方法、装置、终端和存储介质
CN108984851B (zh) 一种带时延估计的加权高斯模型软测量建模方法
CN114386580A (zh) 决策模型训练、决策方法、装置、电子设备及存储介质
CN111191113A (zh) 一种基于边缘计算环境的数据资源需求预测和调整方法
CN116502774B (zh) 一种基于时间序列分解和勒让德投影的时间序列预测方法
CN117407771A (zh) 基于数字孪生的轴承健康状态评估方法、装置及相关设备
Liu et al. Signal denoising method combined with variational mode decomposition, machine learning online optimization and the interval thresholding technique
US20230386448A1 (en) Method of training speech recognition model, electronic device and storage medium
CN116826739A (zh) 基于单变量时序的超短期农业电力负荷预测方法及装置
CN115545168A (zh) 基于注意力机制和循环神经网络的动态QoS预测方法及系统
Chu et al. A new parametric adaptive nonstationarity detector and application
CN112836381B (zh) 一种基于多源信息的船舶剩余寿命预测方法及系统
CN115630561A (zh) 一种神经网络热负荷预测模型自动优化方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant