基于Matlab曲线拟合求取地表沉陷预计参数的程序实现与优化

2018-03-20高超

采矿与岩层控制工程学报 2018年1期

高超

(天地科技股份有限公司开采设计事业部，北京 100013)

煤矿大规模开采造成地表沉陷与裂缝，对生态环境造成较大影响，其影响程度、影响范围及地表沉陷预计是煤矿开采沉陷的主要研究内容之一[1-2]，我国开采沉陷预计方法主要有概率积分、典型曲线法与剖面函数法[3-5]。“三下”采煤沉陷预计无论采用哪种方法，沉陷预计参数的准确性与科学性选取直接关系到地表沉陷预计结果的可靠性，基于地表移动观测站实测资料进行拟合求取概率积分法参数方法仍是最可靠、最科学的参数选取基础来源，并且概率积分法在我国煤矿区应用也最为广泛。

Matlab主要用于算法开发与计算、数据可视化、数据分析、图像处理等领域[6-7]，有一个便于用户使用的视窗环境，并在一定程度上摆脱了非交互式程序设计语言的编辑模式。Matlab内置600多个数学算法函数，并同时支持.m文件应用程序和指令控制；它的指令表达符合数学、工程的思考方式，并能直接调用多种软件程序命令，实现工程应用开发与软件对接；图形处理方面可以进行二维和三维的可视化、图象处理、动画和表达式作图。

开采沉陷学者对于概率积分法曲线拟合求参的方法实现与优化已经进行了大量的探索。吴侃等[8]基于实测资料对概率积分法预计模型进行了修正，解决了预计结果与实测资料在某些方面的不符现象。查剑锋等[9]在综合分析众多实测的基础上，分析了概率积分法地表沉陷预计的误差来源。刘钟森等[10]提出了在概率积分法中引入坐标修正参数，减少了求取预计参数过程中的计算边界坐标假定问题；王云广[11]使用Matlab提供的nlinfit函数实现实测数据的曲线拟合，但nlinfit函数较简单、不如lsqcurvefit函数全面，且对于限定区间的拟合求参精确解有待提高；汪桂生[12]应用C#程序及Matlab软件GUI建立了用户交互界面，实现了矩形单工作面开采地表沉陷实测数据拟合求取程序；谷拴成等[13]对概率积分法进行了适当修正并基于Matlab软件编写了山区浅埋煤层地表沉陷预计参数的求取。

1 基于Matlab的概率积分法曲线拟合求参原理

1.1 应用Matlab曲线拟合求参原理

由地表移动观测站实测资料可获得地表水平移动和下沉数据，并根据不同的求参方法进行求取地表沉陷计算参数，主要有以下3种方法：

(1)根据走向剖面主断面下沉数据求得相关参数，并将一部分参数固定，然后根据该剖面主断面对应的水平移动数据求得其余地表沉陷计算参数。

(2)根据两主断面下沉和水平移动数据使得[VWVW]+[VUVU]=min，一次性求得所有地表沉陷计算参数。

(3)更理想的处理方法是先根据下沉数据拟合求得一组沉陷计算参数，然后再由水平移动数据获得另一组参数，对比两组参数，最后确定出地表沉陷计算参数的综合值。

本文地表移动观测站实测数据拟合求取地表沉陷计算参数即采用了第3种方法。曲线拟合求参的基本原理可表达为：若关于自变量X和待定参数P的函数为Y=f(X，P)，则根据给出(x，y)的n对观测数据(xi，yi)，i=1，2，3，…n，计算出当Q取得最小数值的一组参数P，如式(1)所示：

(1)

拟合过程中利用最小二乘法原理使VTV取最小值，从而选出一组参数P使得拟合曲线与实测结果偏差V的平方和最小。根据下沉数据，以概率积分法为基础用曲线拟合求参法求出一组q，tanβ，S，再利用水平移动数据，求出一组q，b，tanβ，S，最后根据拟合效果并结合现场观测数据完整性、观测效果等情况，确定地表沉陷预计参数的综合值。

曲线拟合求参过程中可使用Matlab软件的非线性最小二乘拟合函数，即nlinfit，lsqcurvefit，lsqnonlin函数，为节约计算时间，这3种拟合函数可给定初始参数值，然后编写概率积分法地表移动变形函数式，在编写好的拟合函数式的基础上计算，计算完成后返回待求的参数组。

1.2 求参过程中概率积分法函数式改进

我国煤矿采掘工程平面图中原始坐标往往是7位或8位数，如果拾取图中的地表移动观测站各测点实测坐标，并将该原始坐标直接应用于概率积分法函数中，对于调用Matlab内置非线性拟合内置函数应用实测或试验数据拟合求取概率积分法计算参数过程中，往往采用给定参数区间、给定初始值、设定计算精度的迭代法进行计算；当原始坐标位数较多时，计算内存占用量较大，函数的执行效率将降低。

当对实测数据或实验数据拟合求取地表移动预计参数完成后，将主要影响半径r值再乘以1000，继而应用该值进行后续的预计参数求取。

2 Matlab曲线拟合求参程序编写

2.1 下沉数据曲线拟合求参

试验矿井位于内蒙古东胜煤田北部，试验工作面地表地形属高原侵蚀性丘陵地貌，大部分地区为低矮山丘，第四系松散层分布较为广泛，植被稀疏，为半荒漠地区。试验工作面四邻均为未开采区域，工作面煤层平均厚度3.4m，倾角1～5°，平均埋藏深度120m；工作面开采长度约1782m，工作面斜长270m，采用走向长壁式采煤方法，全部垮落法管理顶板。根据概率积分法函数公式，应用Matlab程序编写由实测数据的概率积分法预计参数拟合求取程序。

2.1.1 倾向主断面求参程序编写

利用Matlab中的xlsread直接读取指定存储目录下Excel表格中倾向全盆地观测线的观测数据，应用如下命令：

X=xlsread('C:MATLAB7work est2高程.xlsx'，'sheet1'，'d2:d42')；

同理应用该命令读取各次观测的高程值；

Y3=1000*(Y2-Y1);%实验或实测下沉值转换

subs=inline('-3400/2*beta(1)*(erf((X/1000-beta(3)-0).*(sqrt(pi)/beta(2)))+1)+3400/2*beta(1)*(erf(((X/1000-(275+beta(4))+0)*sin((beta(5)+1)*pi/180)/sin((beta(5))*pi/180)).*(sqrt(pi)/beta(2)))+1)'，'beta'，'X');%编写下沉预计函数

parm=lsqcurvefit(subs，[0.7，50，5，5，50]，X，Y3，，，opt);%赋参数的初始值

H=120;%煤层平均埋藏深度

beta(1)=parm(1);beta(2)=parm(2);beta(3)=parm(3);beta(4)=parm(4);beta(5)=parm(5);%拟合结果赋值给概率积分公式参数

q= beta(1);tanb=H/(beta(2)*1000);S1=beta(3);S2=parm(4);the=parm(5);%将拟合参数赋值到沉陷预计参数表示

parameters1=[q，tanb1，S3，S4 ，tanb2] %——输出预计参数

W=-3400/2*beta(1)*(erf((X/1000-beta(3)-0).*(sqrt(pi)/beta(2)))+1)+3400/2*beta(1)*(erf(((X/1000-(275+beta(4))+0)*sin((beta(5)+1)*pi/180)/sin((beta(5))*pi/180)).*(sqrt(pi)/beta(2)))+1); %将输出参数代入预计函数

plot(X，Y3，'o'，X，W，'r')；%绘制下沉曲线

xlabel('倾向L (m)');ylabel('下沉W (mm)'); %坐标轴输出控制

set(gca，'XTick'，-250:50:450);set(gca，'YTick'，-3000:500:500)；%坐标轴刻度控制

legend('原始数据'，'拟合曲线');%图例及图形设置

2.1.2 走向主断面变采高求参程序及改进

根据试验矿井生产调度日报表及井上下对照图，沿试验工作面推进方向的观测点下沉主要在2012年8—10月，其中7月份平均采高3.25m，8月份平均采高3.25m，9月份平均采高3.40m，10月份平均采高3.53m，11月份平均采高3.53m。

煤矿开采过程中的各个月份的平均采高并非保持不变，而是变化、浮动的，对于沿工作面推进方向的地表移动观测站实测数据曲线拟合求参，通常将工作面的平均采高代入地表沉陷预计公式中，也有部分学者采用垂直划分、将采掘工作面划分为若干个月份的采掘之和，这对于采高变化大的工作面拟合求参来说误差较大。为解决沿工作面推进方向不同位置采厚的变化，作者提出将实际采掘工作面以水平面划分为多个不同煤层厚度的工作面叠加的原理来进行地表沉陷预计参数求取的方法，减小了工作面平均采高法与垂直划分法的预计参数求取误差，尤其解决了工作面垂直划分法求取参数中拐点偏移距的离散性问题。

该试验工作面7—11月份最小采高为3.25m、最大采高为3.53m，并将工作面沿水平面划分为3层：最下层的3.25m煤层、中间层的0.15m煤层与最上层的0.13m煤层，依据工作面实际回采进度，沿工作面推进方向变采高的下沉预计函数(不考虑坐标位数对计算速度影响)可表示为：

W=inline('-3250/2*beta(1)*(erf((X-beta(3)).*(sqrt(pi)/beta(2)))+1)+3250/2*beta(1)*(erf((X-1780+beta(4)).*(sqrt(pi)/beta(5)))+1)-150/2*beta(1)*(erf((X-233-beta(3)).*(sqrt(pi)/beta(2)))+1)+150/2*beta(1)*(erf((X-233-1780+beta(4)).*(sqrt(pi)/beta(5)))+1)-130/2*beta(1)*(erf((X-471-beta(3)).*(sqrt(pi)/beta(2)))+1)+130/2*beta(1)*(erf((X-471-1780+beta(4)).*(sqrt(pi)/beta(5)))+1)'，'beta'，'X');%编写走向盆地观测线下沉预计函数

2.2 水平移动数据曲线拟合求参

同理，根据概率积分法函数式，应用Matlab程序编写实测数据的概率积分法拟合求参程序：

disp=inline('2448*beta(1)*(exp(-pi*((X/1000-beta(3)).^2)./(beta(2).^2))-exp(-pi*((X/1000-(270-beta(4))).^2)./(beta(2)^2)))-2448/2*(erf((X/1000-beta(3)).*(sqrt(pi)/beta(2)))-erf((X/1000-(270-beta(4))).*(sqrt(pi)/beta(2))))*cos(beta(5)*pi/180)/sin(beta(5)*pi/180)'，'beta'，'X')；%水平移动预计函数

parm=lsqcurvefit(disp，[0.3，60，5，5]，X，Y，[0.20，30，-15，-15]，[0.55，90，15，15]，opt)；%赋初始值

parameters1=[b，tanb，S3，S4]；%输出预计参数

Disp=2448*beta(1)*(exp(-pi*((X/1000-beta(3)).^2)./(beta(2).^2))-exp(-pi*((X/1000-(270-beta(4))).^2)./(beta(2)^2)))-2448/2*(erf((X/1000-beta(3)).*(sqrt(pi)/beta(2)))-erf((X/1000-(270-beta(4))).*(sqrt(pi)/beta(2))))*cos(beta(5)*pi/180)/sin(beta(5)*pi/180);%将输出参数代入预计函数

plot(X，Y3，'o'，X，Disp，'r')；%绘制水平移动曲线

2.3 曲线拟合精度检验

根据测量平差及误差理论对实测的高程和平面观测数据用给定拟合曲线模型进行求取地表沉陷计算参数时，往往利用最小二乘法使偏差平方和最小，在根据该原理进行拟合效果的判定时，有时拟合效果的好坏不能很好地以该原理反映出来，因为测点相对于整条测线范围只是测线上的部分数据，应归类为样本数据，为了更好地进行曲线拟合效果的判别，将测点残差平方和Q和样本点数结合起来，引入样本均方差进行曲线拟合效果的检验。样本均方差S的定义如式(2)所示，式中N为实测点个数。

(2)

对比分析各个测线的概率积分法拟合曲线及样本均方差可知，各测线样本均方差较小，即拟合效果较好，样本均方差较大则拟合效果一般。

3 地表沉陷计算参数求取应用

基于Matlab平台编写的曲线拟合程序对试验煤矿地表观测站实测下沉与水平移动数据进行曲线拟合求参。曲线拟合求参过程中，沿煤层走向方向下沉曲线和水平移动曲线，将坐标原点定在开切眼侧的采空区边界处(X轴沿工作面推进方向为正)；沿煤层倾斜方向的下沉曲线和水平移动曲线，将坐标原点定在下山侧的采空区边界处(X轴沿煤层上山方向为正)；拐点偏移距若为正值，表示在采空区一侧，反之表示在煤柱一侧。

以试验工作面实测地表移动变形数据拟合求参为例，倾向测线实测下沉数据对地表移动预计参数拟合结果如图1所示，应用倾向测线实测水平移动数据对地表移动预计参数拟合结果如图2所示，应用走向测线变采高实测下沉数据对地表移动预计参数拟合结果如图3所示。

图1 倾向测线下沉数据拟合求参曲线

图2 倾向测线水平移动数据拟合求参曲线

图3 走向测线变采高下沉数据拟合求参曲线

将由实测数据应用概率积分法拟合求得的结果与样本标准差整理进入表1中，并最终得出了该试验矿井地质采矿条件下的地表沉陷预计参数。

试验工作面地表移动观测站倾向观测线实测最大下沉值为2448mm、最大水平移动值为430mm，走向观测线的实测测点的最大下沉值为2571mm。由拟合结果可知倾向测线下沉数据拟合求参样本标准差为26.2mm，约为最大下沉值的1.1%，说明该测线的下沉数据很符合概率积分法预计公式，拟合求取的概率积分法参数效果很好；倾向测线水平移动数据拟合求参样本标准差为71.6mm，约为最大水平移动值的16.6%(倾向测线初次平面观测时间稍滞后，致使水平移动数据拟合求参效果较差)，说明该测线的水平移动数据比较符合概率积分法预计公式，拟合求取的概率积分法参数效果一般；走向测线下沉数据拟合求参样本标准差为29.9mm，约为最大下沉值的1.2%，说明对于变采高工作面将工作面沿水平面划分并拟合求取概率积分法计算参数的方法与结果很好。