CN103258112B - 一种基于mde和分形优化svm的周期来压预测方法 - Google Patents

Abstract

本发明公开了一种基于MDE和分形优化SVM的周期来压预测方法，<b />主要步骤包括对已知支架周期来压荷载曲线使用多重差异进化算法（MDE）进行拟合，将每重拟合形成的单一正弦曲线与上次差余曲线（）再作差余曲线()。将这些图通过分形几何的盒子法计算维度和相关系数（r）。将每条的维度、<i>r</i>和支架相对距离(L)作为输入值，对应的的周期、缩放系数和纵移系数作为目标值，使用支持向量机（SVM）进行训练。通过对维度和r规律的研究得到拟设置支架处荷载各的维度和r，带入训练后的SVM模拟得到的、和，进而得到的表达式。将上述求和即为所求拟设置支架处的周期来压荷载。<b />可广泛用于液压支架周期来压拟合和预测。

Description

一种基于MDE和分形优化SVM的周期来压预测方法

技术领域

本发明涉及地下采矿工程周期荷载预测问题，特别是涉及基于MDE和分形优化SVM的周期来压预测方法。

背景技术

在巷道内的液压支架是承受顶板压力的主要构件，在设置支架时要考虑很多因素，其中荷载的形式和变化规律是主要的考察因素。如何了解将要设置液压支架处的周期来压荷载形式和变化规律成为关键问题。当然，这个问题在未设置液压支架前是不可测的，这样只能通过预测来解决上述问题。

现有技术中没有解决在不可测情况下，如何预测将要设置液压支架所受到的荷载形式和变化规律的问题。通过已知的液压支架的荷载表示未设置液压支架的荷载，由于未设置液压支架处的荷载不可测，无法进行分析。在这种情况下，根据已有的研究和模型的分析，可以确定以设置的液压支架之间的荷载形式和变化规律存在着一些联系，而且他们的荷载是可测的。根据模型对不同相对位置支架周期来压荷载图像的分析，发现就支架相对位置而言，差余曲线图形的维数和相关系数r有一定的线性关系。根据所求支架与已知支架的相对位置关系和该线性关系即可求得所求支架周期来压荷载图像各对应差余曲线的维数和相关系数r。这样就可将所求支架荷载与已知支架荷载联系起来，再通过一定方法进行映射。

模型通过MDE得到周期来压荷载的正弦波分解形式；分形理论计算了各差余曲线的维数和相关系数r；通过SVM模拟出未设置液压支架处荷载曲线各正弦项的参数；将各项叠加得到该荷载的预测表达式。通过与实际监测的结果对比表明了结果的正确性及方法的正确性，本发明涉及的周期来压预测方法能够解为顶板支护提供合理依据。

发明内容

1.多重差异进化算法拟合

设原始测量的液压支架荷载数据的曲线为Y^O；对Y^O使用差异进化算法拟合得到的正弦曲线集合为y_i＝S_i×sin(2πx/T_i-π/2)+D_i，S_i,T_i,D_i是拟合求得的参数(用TSM，TZM，TDM分别表示S_i,T_i,D_i各自的矩阵)，分别表示周期(T_i)、缩放系数(S_i)和纵移系数(D_i)；i∈[1,n]，n为计算的代数，即拟合的次数(拟合曲线项数)。定义差余曲线E_i为：本代差余曲线是上代差余曲线与本代数据反演得到的正弦曲线y_i的差的曲线，如E₁＝Y^O、E₂＝Y^O-y₁、E₃＝E₂-y₁…。那么对Y^O的拟合最终结果就是并通过MATLAB实现。

使用差异进化算法有两个目的，第一个是反演出S_i,T_i,D_i作为SVM训练时的目标值；另一个是得到y_i，进而得到E_i的图像，作为分形几何计算维度和r的基础。

2.盒子法计算E_i图像的维数和相关系数r值

盒子法是利用覆盖法测量分维曲线的最基本方法之一，实现简单。通过对标准曲线的测定，结果与理论计算值一致性良好，较适合平面曲线的分维测量。该方法用边长为1的方盒覆盖曲线,将该方盒分割成含有2ⁿ个小方盒的网格集，小方盒的边长为2^-n，用这个网格集覆盖轮廓曲线，统计出含有轮廓片段的小方盒数量为M(n)，则曲线的分形维数为:

$F=\underset{n\&RightArrow;\mathrm{\&infin;}}{\lim }\frac{\lg M\left(n\right)}{n\lg 2}$

鉴于处理二维折线图像，模型使用盒子法计算E_i图像的维数和r值。使用FractalFox2.0对形成的E_i图像进行处理。维数、r和L作为SVM训练的输入值。

3.支持向量机进行预测

使用Epsilon回归算法来进行训练和模拟。前文提到，维数、r和L作为训练的输入值，S_i,T_i,D_i作为目标值，进行训练。模拟的输入值是所求支架荷载的维数、r和L，输出值为所求荷载曲线的正弦项y_i的S_i,T_i,D_i。对于没有数据的所求支架荷载的E_i图像维数和r是通过线性拟合实现的，详见实例的具体计算过程。

4.模型算法流程

上文已经介绍了相关的方法以及在模型算法中的作用。下面用流程图的形式表示算法，如图1所示。

附图说明

图1方法流程图

注：长方形表示数据处理动作，斜四边形表示原始数据，箭头表示数据及其传递方向。

图210#、35#、60#、85#的85次循环左右两支架的荷载变化图

图38条曲线与其拟合曲线的对照图

图4Y^10#z的第2带差余曲线图图5图4对应的分形特征曲线图

注：Y^10#z的差余曲线一共20个图，即这里只列出了图。

图5图4对应的分形特征曲线图

注：Y^10#z的差余曲线一共20个图，即这里只列出了图。

图68条拟合后曲线各项参数图

图78条曲线第一代差余曲线维度与相对距离(L)的关系

图8107#支架模拟荷载与实际荷载对比图

具体实施方式

为使本发明的上述目的、特征和优点更加明显易懂，下面结合使用到的相关理论和具体实施方式对本发明作进一步详细的说明。

实施例为某矿二采区1212回采工作面的周期来压，其位于北二12煤集运巷以西，北邻北二采区1210-1采空区，南侧为未采区。工作面设计采高4.0米，沿顶板推进。采用倾斜长壁后退式综合机械化采煤方法。进刀方式为端部斜切进刀，返往一次割两刀煤，循环进尺为0.8m。

1212工作面开切规格为8.4m×3.5m,采用8m锚索、金属网、钢带联合支护。工作面开采初期采用单向割煤，从尾往头推溜。一是为了调整运输机与转载机合理搭接长度。二是为了保证支架平行进入煤壁，支架顶梁进入煤壁后开始对支架大柱压力表进行读数记录。

在工作面建五条测线监测支架工作阻力变化情况，测线位置分别布置在10#、35#、60#、85#、107#左右支架，每个循环未对支架柱压力表读数一次，并作为记录结果。

原始数据如图2所示，为10#、35#、60#、85#的85次循环左右两支架的荷载变化图。

对图2中8条曲线进行拟合，拟合的方法是MDE。8条实际数据曲线分别为Y^10#z、Y^10#y、Y^35#z、Y^35#y、Y^60#z、Y^60#y、Y^85#z、Y^85#y，其中z表示左，y表示右，y^10#z、y^10#y、y^35#z、y^35#y、y^60#z、y^60#y、y^85#z、y^85#y表示对应的拟合曲线。现以Y^10#z为例说明拟合函数y^10#z的建立过程。基于傅里叶拟合思想和8条曲线都表现为某种正弦曲线的形式，将Y^10#z分解为多项正弦函数的和，如式(1)所示。

$y^{10\#z}=\underset{i=n}{\overset{1}{\&Sigma;}}{\mathrm{TSM}}_{10\#z,i}\&times;\sin (\frac{2\mathrm{\&pi;}x}{{\mathrm{TZM}}_{10\#z,i}}-\frac{\mathrm{\&pi;}}{2})+{\mathrm{TDM}}_{10\#z,i}---\left(1\right)$

式中：TZM_10#z,i表示周期，对y^10#z第i代拟合(第i项)的正弦曲线的周期；TSM_10#z,i表示缩放系数，对y^10#z第i代拟合(第i项)的正弦曲线的振幅程度；TDM_10#z,i表示纵移系数，对y^10#z第i代拟合(第i项)的正弦曲线的纵向平移量；x为掘进深度/m。TZM，TSM，TDM分别表示周期T_i、缩放系数S_i和纵移系数D_i矩阵。

使用差异进化(DE)进行反演的方法，通过计算y^10#z的第一项进行说明，第一项如式(2)所示。

$y_1^{10\#z}={\mathrm{TSM}}_{10\#z,1}\&times;\sin (\frac{2\mathrm{\&pi;}x}{{\mathrm{TZM}}_{10\#z,i}}-\frac{\mathrm{\&pi;}}{2})+{\mathrm{TDM}}_{10\#z,1}---\left(2\right)$

参数识别过程如下:

1)收集现场检测数据Y^10#z，进行数据处理。对数据有效性进行分析，形成适合的数据结构。

2)选择正弦曲线作为反演分析数学模型。

3)Y^10#z、x作为回归函数拟合数据，选择收敛最快且值最小时的缩放因子F和交叉概率Cr。

4)利用确定好的缩放因子F和交叉概率Cr和实际数据，建立二维反演参数识别与优化模型。

5)以参数TSM_10#z,1、TZM_10#z,1和TDM_10#z,1作为优化变量，构建回归模型与现场监测数据Y^10#z的最小二乘函数关系进行DE运算，即设定DE初值，按照DE步骤进行参数搜索。

6)算法达到最佳收敛后，即可输出获得的TSM_10#z,1、TZM_10#z,1和TDM_10#z,1参数。

要注意的是计算第二项时第3步和第5步的Y^10#z改为以此类推。这样就可以求得y^10#z，进而可以求得8条拟合曲线。关于拟合曲线的项数(差余曲线代数)还确定不了，先计算20项，然后观察其收敛值(msumM)情况。如表1所示为y^10#z的参数；8条曲线拟合后如图3所示。拟合质量如表2所示。

表1y^10#z的参数表

注：msumM表示差余曲线上各点值的平方和。

表28条曲线的拟合质量

注：每条曲线计算20项，但是有的曲线提前收敛，或到某项后收敛值变大，这时在它们后面的项舍掉。

图3显示了曲线的对比情况，从图中看出虽然有些部分的曲线符合度不佳，但是曲线的走势基本一致，拟合曲线能按照原曲线的变化而变化。10#支架左、10#支架右和60#支架右的对比图质量较差，主要原因是由于拟合的项数较小，无法体现更丰富的变化，但是其msumM最小项后的项以失去意义应舍掉。总体拟合质量如表2所示，拟合质量可以达到要求，而且这里主要关注的是变化趋势和程度，因为这些决定了图像的分形维数和r值。

上述过程得到了两个方面的数据，一方面是8条拟合曲线的TSM、TZM和TDM矩阵，另一方面是E_i差余曲线。

E_i是为了使用分形几何计算参数和r设置的，为了解决已设与未设支架处荷载不能建立数值关系的缺点，而借助图像分析建立分形几何关系。就Y^10#z曲线的E_i进行说明，第一代差余曲线第二代用于分析的Y^10#z差余曲线图为如图4所示为

使用FractalFox2.0对20×8条曲线进行处理得到图的维数和r，参数为处理方式：盒子法(BoxCounting)、FromBoxSize:2、ToBoxSize:100、StepSize:1。计算得到各曲线的各E_i的维数和r。y^10#z的参数如表1所示，分形特征曲线如图5所示。

专利解决问题的关键是建立在没有数据情况下，依靠图像变化规律找到E_i图像内部和图像之间的联系。为了进行下一步的支持向量机的训练和模拟，将8条曲线得到的参数值归一到[-1,1]，如图6所示为归一后y^10#z各项参数变化图。

表1中可以看出，随着msumM的减小，即曲线拟合质量增加时曲线分形几何维度增加。这说明曲线变得粗糙了，但是更接近于y＝0的平衡位置，r随着维数的变化也有相应的变化。图6中，基本上每个图内的各参数曲线之间是有规律的，图之间相同参数曲线变化基本相同。说明在支架荷载的变化曲线中，普遍存在着某种关系，这个关系可以通过该方法计算得到的维数和r进行描述。另一方面，已知支架和未设立的支架存在着相对位置关系，按照支架之间间隔1m计算，10#、35#、60#、85#、107#相对于107#的距离分别为97m、72m、47m、22m、0m。随着距离的不同拟合曲线对应各项E_i的维数显现出线性规律，如图7所示为左右两侧的四条第一代差余曲线E₁维度。

图7的维数和相对距离的线性关系比较明显，后代差余曲线维数也有形同的规律。对于另一个分形参数相关系数r，其与msumM虽然没有完全符合的对应关系，但是其变化与维数的变化频率和方向基本相同，也可以看做是与msumM的对应，且与相对距离也存在像维数那样的关系。综上，维数和r可以反映不同相对距离支架荷载曲线之间的关系，且自身的变化与相对距离也有关系，所以维数、r和L就是已知与未知荷载曲线之间的桥梁。

支持向量机预测首先要有三个条件：

1)训练的输入值。上文提到的(L、维数、r)_i矩阵即为输入值。首先要确定预测的107#支架荷载曲线的项数(代数)，该项数与10#、35#、60#、85#的拟合曲线的项数相对应。表2中列出了各曲线的项数，最小为6，最大为20。如果选择6项，预测后的叠加结果的变化就比较简单，曲线拟合效果较差；如果选择20，那么有4个拟合曲线的有效项不足20，扩展到20项拟合曲线就不准确。综合考虑，选择15项效果较好。那么8个拟合曲线的项数都取15项。对于不到15项的拟合函数，通过添加正反项的方式补充到15项，如y^10#y只有6项，那么缺9项。这9项的所有参数与第6项相同，只是第7项的TSM与第6项相同，第8项的TSM与第6项相反，依次类推。这样做可以最大程度的相互抵消多余项，使msumM变化最小，保证函数图像的稳定性。

模拟和预测都是左线和右线支架分开进行的，因为其函数图像的特征不同，各自单侧才与相对距离有对应关系。所以(L、维数、r)_i ^Z为左线训练输入值，3×60矩阵；(L、维数、r)_i ^Y为右线训练输入值，3×60矩阵；

2)训练目标值。因为目标值周期，缩放系数和纵移系数与训练数据构成的支持向量机结构是不同的，所以将目标值定义为(周期)_i ^Z、(缩放系数)_i ^Z、(纵移系数)_i ^Z、(周期)_i ^Y、(缩放系数)_i ^Y、(纵移系数)_i ^Y，6个1×60的矩阵。

3)模拟输入值。模拟输入值为δ(L、维数、r)_i。这些只是对于107#支架而言的，所以相对距离r＝0。由图7及分析可知，可通过线性拟合得到以10#、35#、60#、85#数据为基础的107#支架荷载拟合曲线各项的维数和r。使用MATLAB的polyfit函数拟合并计算。得到的结果如表3所示。

表3拟合得到的107#左右支架的维数和r

那么δ(L、维数、r)_i ^Z，3×15矩阵，为模拟左线的输入值；δ(L、维数、r)_i ^Y，3×15矩阵，为模拟右线的输入值。

当然上述SVM相关数据都要进行归一处理才能满足使用要求。将上述值带入SVM进行训练和模拟，参数设置如下：使用Epsilon回归算法；拉格朗日乘子上界C＝100；不敏感损失函数的参数e＝0.2；ker＝struct('type','gauss','width',0.6)。得到模拟输出值后进行逆归一运算，得到真实值如表4所示。

表4模拟得到的107#左右两线荷载拟合曲线的各项参数

根据式(1)得到式(3)和(4)：

$y^{107\#z}=\underset{i=n}{\overset{1}{\&Sigma;}}{\mathrm{TSM}}_{107\#z,i}\&times;\sin (\frac{2\mathrm{\&pi;}x}{{\mathrm{TZM}}_{107\#z,i}}-\frac{\mathrm{\&pi;}}{2})+{\mathrm{TDM}}_{107\#z,i},i=\&lsqb;1,\mathrm{...},15\&rsqb;---\left(3\right)$

$y^{107\#y}=\underset{i=n}{\overset{1}{\&Sigma;}}{\mathrm{TSM}}_{107\#y,i}\&times;\sin (\frac{2\mathrm{\&pi;}x}{{\mathrm{TZM}}_{107\#y,i}}-\frac{\mathrm{\&pi;}}{2})+{\mathrm{TDM}}_{107\#y,i},i=\&lsqb;1,\mathrm{...},15\&rsqb;---\left(4\right)$

将表4中数据带入得到的模拟曲线并与实际107#左右支架荷载曲线对比如图8所示。

图8中左右两图的均方差和最大误差分别是8.8267、1.3504和6.2610、0.7982。模拟质量基本达到要求。更为关键的是在图8中，可以明显的看到模拟曲线跟随着实际监测数据的变化而变化，这对预测支架周期来压荷载变化有重要意义。同时从振幅上看，模拟曲线基本包括了实际荷载曲线，即其振幅值可以大体反映来压的最大值。综上两点可以说明，虽然模拟质量一般，但是本方法可以反映支架荷载的基本形式和变化规律。特别是在未开挖及未设置支架处，无法对数据进行测量，但已知相邻多个支架荷载变化的情况下，该方法显得尤为有效。

Claims (4)

1.一种基于MDE和分形优化SVM的周期来压预测方法，主要步骤包括：

步骤一：对已知支架周期来压荷载曲线使用多重差异进化算法(MDE)进行拟合，将每重拟合形成的单一正弦曲线与上次差余曲线E_i再作差余曲线E_i+1，定义差余曲线E_i为：本代差余曲线是上代差余曲线与本代数据反演得到的正弦曲线y_i的差的曲线；步骤二：将这些E_i图通过分形几何的盒子法计算维度和相关系数r；

步骤三：将每条E_i的维度、相关系数r和支架相对距离L作为输入值，对应的E_i的周期T_i、缩放系数S_i和纵移系数D_i作为目标值，使用支持向量机(SVM)进行训练；

步骤四：通过对维度和相关系数r规律的研究得到拟设置支架处荷载各E_i的维度和相关系数r，带入训练后的SVM模拟得到y_i的T_i、S_i和D_i，进而得到y_i的表达式；设原始测量的液压支架荷载数据的曲线为Y^O；对Y^O使用差异进化算法拟合得到的正弦曲线集合为y_i＝S_i×sin(2πx/T_i-π/2)+D_i，S_i,T_i,D_i是拟合求得的参数；S_i,T_i,D_i分别表示E_i的周期、缩放系数和纵移系数，i∈[1,n]，n为计算的代数，即拟合的次数；

步骤五：将上述y_i求和即为所求拟设置支架处的周期来压荷载。

2.根据权利要求1所述一种基于MDE和分形优化SVM的周期来压预测方法，其特征在于，10#左线支架的荷载正弦曲线y^10#z分解为多项正弦函数的和： $y^{10\#z}=\underset{i=n}{\overset{1}{\&Sigma;}}{\mathrm{TSM}}_{10\#z,i}\&times;sin(\frac{2\mathrm{\&pi;}x}{{\mathrm{TZM}}_{10\#z,i}}-\frac{\mathrm{\&pi;}}{2})+{\mathrm{TDM}}_{10\#z,i}$ 式中：TZM_10#z,i表示周期，对y^10#z第i代拟合的正弦曲线的周期；TSM_10#z,i表示缩放系数，对y^10#z第i代拟合的正弦曲线的振幅程度；TDM_10#z,i表示纵移系数，对y^10#z第i代拟合的正弦曲线的纵向平移量；x为掘进深度，单位为m；TZM，TSM，TDM分别表示周期T_i、缩放系数S_i和纵移系数D_i的矩阵；使用差异进化进行反演的方法，y^10#z的第一项计算如下： $y_1^{10\#z}={\mathrm{TSM}}_{10\#z,1}\&times;\sin (\frac{2\mathrm{\&pi;}x}{{\mathrm{TZM}}_{10\#z,1}}-\frac{\mathrm{\&pi;}}{2})+{\mathrm{TDM}}_{10\#z,1}$ 参数识别过程如下：1)收集现场检测数据y^10#z，进行数据处理；对数据有效性进行分析，形成适合的数据结构；2)选择正弦曲线作为反演分析数学模型；3)y^10#z、x作为回归函数拟合数据，选择收敛最快且值最小时的缩放因子F和交叉概率Cr；4)利用确定好的F和Cr和实际数据，建立二维反演参数识别与优化模型；5)以参数TSM_10#z,1、TZM_10#z,1和TDM_10#z,1作为优化变量，构建回归模型与现场监测数据y^10#z的最小二乘函数关系进行差异进化运算，即设定差异进化初值，按照差异进化步骤进行参数搜索；6)算法达到最佳收敛后，即可输出获得的TSM_10#z,1、TZM_10#z,1和TDM_10#z,1参数；其中，10#支架代表序号为10的支架；以此类推，计算各项的y_i ^10#z，可以求得y^10#z。

3.根据权利要求1所述的一种基于MDE和分形优化SVM的周期来压预测方法，其特征在于，所述差余曲线的构造，分形几何计算参数和相关系数r的设置，使用FractalFox2.0对20×8条曲线进行处理得到图的维数和r，参数处理方式：盒子法(BoxCounting)、FromBoxSize:2、ToBoxSize:100、StepSize:1，计算得到各曲线的各E_i的维数和r；SVM模拟使用Epsilon回归算法来进行训练和模拟，维数、r和L作为训练的输入值，S_i,T_i,D_i作为目标值，进行训练；模拟的输入值是所求支架荷载的维数、r和L，输出值为所求荷载曲线的正弦项y_i的S_i,T_i,D_i；对于没有数据的所求支架荷载的E_i图像维数和r是通过线性拟合实现的；使用Epsilon回归算法，其中，拉格朗日乘子上界C＝100；不敏感损失函数的参数e＝0.2；编程结构ker＝struct('type','gauss','width',0.6)。

4.根据权利要求1所述的一种基于MDE和分形优化SVM的周期来压预测方法，其特征在于，要输入训练的输入值，以107#支架为实际处理对象，首先要确定预测的107#支架荷载曲线的项数，该项数与10#、35#、60#、85#的拟合曲线的项数相对应，项数确定的方法如下：各曲线的项数，最小为6，最大为20；选择6项，预测后的叠加结果的变化简单，曲线拟合效果差；选择20，那么有4个拟合曲线的有效项不足20，扩展到20项拟合曲线就不准确；综合考虑，选择15项，那么8个拟合曲线的项数都取15项，对于不到15项的拟合函数，通过添加正反项的方式补充到15项；模拟和预测都是左线和右线支架分开进行的，因为其函数图像的特征不同，各自单侧才与相对距离有对应关系，所以(L、维数、r)_i ^Z为左线训练输入值，为3×60矩阵；(L、维数、r)_i ^Y为右线训练输入值，为3×60矩阵；所述的训练的目标值定义为(周期)_i ^Z、(缩放系数)_i ^Z、(纵移系数)_i ^Z、(周期)_i ^Y、(缩放系数)_i ^Y、(纵移系数)_i ^Y，6个1×60的矩阵；其中()^Z代表左线的值，()^Y代表右线的值；n#支架代表序号为n的支架。