时间：2013年07月13日 分类：推荐论文 次数：
摘要：本文利用MATLAB 工具箱函数，设计了求取地表移动变形预计参数的程序，结合工作面任意点的观测数据和地质采矿条件，能够很好的求取地表移动变形预计参数。本文通过实例，对所拟合的参数进行检验，证明了此方法的可行性。
关键字：MATLAB,非线性最小二乘法,概率积分法,预计参数
1 引言
MATLAB有美国MathWorks公司推出，其与Mathematica、Maple并称为三大数学软件，具有优秀的数值计算能力和卓越的数据可视化功能。 MATLAB自问世以来，已经发展称为适合多学科、多种工作平台的功能强大的大型软件，广泛应用于大学科研、工程计算等领域，尤其在工程界。
概率积分法通过积分而得到下沉曲线表达式，因其所用的移动和变形预计公式中含有概率积分而得名。这种方法是将矿山岩层移动作为一种服从统计规律的随机现象来讨论。因此，此种方法是以随即介质理论为基础的一种预计方法。概率积分法自20世纪60年代引入我国以来，在许多矿山得到广泛使用，并不断发展，目前是我国进行地表移动变形预计主要的方法之一。地表移动变形预计参数对开采沉陷预计的准确性具有重要意义，因而本文以概率积分法数学模型为基础，借助MATLAB平台，编写求取地表移动变形参数的程序。
2 概率积分法求参的数学模型
利用概率积分法关于任意点的地表下沉预计公式作为求取预计参数的基本数学模型，即：
地表移动观测站任何一个测点的下沉实测值W全能表达成自变量(x,y)(观测点的水平坐标)和概率积分法预计参数q、tgβ、s3、s4、tgβ1、tgβ2、s1、s2、θ的函数，见下式：
假设公式(3)里的参数q、tgβ、s3、s4、tgβ1、tgβ2、s1、s2、ctgθ，已经求出，那么上面的数学公式中就只剩下水平移动系数b这一待求参数。
3 Matlab曲线拟合求参函数
Matlab强大的绘图和用户自定义函数等功能，在数据处理领域被广泛应用，其中当然包括对实验、实测数据的拟合工作。鉴于本文研究内容，着重介绍Matlab中非线性曲线拟合(非线性最小二乘拟合)函数。
Matlab中提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。lsqcurvefit的定义如下：
已知拟合的数据点为：xdata和ydata，lsqcurvefit用以求含参量variable(向量)的向量值函数F(x,xdata)中的参变量x(向量)，使满足下式：
完整调用格式为：
式中：x是拟合的参数最优解，norm是最优解的误差的平方和，res为误差向量，ef是程序结束时的状态指示：ef>0表示结果收敛，ef=0表示函数调用次数或迭代次数达到最大值(该值在options中指定)，ef<0表示不收敛，out则包含数据的一个结构变量，包含实际的迭代次数等，lam是上下界所对应的Lagrange乘子，jac是结果(x点)处的雅可比矩阵，x0为初始参数(缺省时程序自动取x0=0)fun是给出的目标函数，xdata和ydata为拟合的数据，v1,v2是拟合参数的上下界，可以用两个[]代替，options是包含算法控制参数的结构设定(或显示)控制参数的命令为Optimset函数。
Matlab可以很快速的直接读写txt格式和Excel格式的数据文件，所以给开采沉陷数据处理带来了方便，解决了求参和预计过程中繁多的数据手工输入和输出的问题。在Matlab中读写Excel文件可以使用xlsread函数和xlswrite函数，xlsread函数的使用格式为：[num txt]=xlsread(‘filename’，sheet，’range’)，读一个Excel文件名为filename的文件，sheet用来指定页，range用来指定区域，最后把指定的数值读到num中，所有字符串到字符串单元数组txt中。使用xlswrite函数写Excel文件，经常用到的格式为：xlswrite (‘filename’，M，sheet，’range’)，写矩阵或字符串单元数组M到filename中的指定页和指定区域。
4 地表移动变形参数求取与检验
4.1地表移动变形预计参数求取的算法流程
4.2求取沉陷预计参数的程序设计
本文主要用到的两个自定义函数budgets和horizonal_budgets，函数budgets用于求取沉陷预计参数q、tgβ、s3、s4、tgβ1、tgβ2、s1、s2、θ，部分源代码如下：
options=optimset('Display','iter','MaxFunEvals',2^8,'MaxIter',2^14,'TolX',1e-3,'TolFun',1e-4,'LargeScale','off','LevenbergMarquardt','on','LineSearchType','cubicpoly');
f=inline('1./(a(1)*x(3,4)*cos(x(3,5))).*((a(1)*x(3,4)*cos(x(3,5)))/2.*(erf(sqrt(pi).*x(1,:).*a(2)./x(3,1))-erf(sqrt(pi).*(x(1,:)-(x(3,6)-a(5)-a(6))).*a(2)./x(3,1)))).*((a(1)*x(3,4)*cos(x(3,5)))/2.*(erf(sqrt(pi).*x(2,:).*a(3)./x(3,2))-erf(sqrt(pi).*(x(2,:)-(x(3,7)-a(7)-a(8)).*(sin(a(9)+x(3,5))/sin(a(9)))).*a(4)./x(3,3))))','a','x');
[a,norm,res,ef,out,lam,jac]=lsqcurvefit(f,a,[x;y;z],w,[],[],options);
函数horizonal_budgets用于求取沉陷预计参数b，部分源代码如下：
options=optimset('Display','iter','MaxFunEvals',2^8,'MaxIter',2^14,'TolX',1e-8,'TolFun',1e-12,'LargeScale','on','LevenbergMarquardt','on')
f=inline('b.*(x(4,1).*x(3,4).*cos(x(3,5))).*(exp(-pi.*(x(1,:).*x(4,2)./x(3,1)).^2)-exp(-pi.*((x(1,:)-(x(3,6)-x(4,5)-x(4,6))).*x(4,2)./x(3,1)).^2)).*(x(4,1)*x(3,4)*cos(x(3,5))/2.*(erf(sqrt(pi).*x(2,:).*x(4,3)./x(3,2))-erf(sqrt(pi).*(x(2,:)-(x(3,7)-x(4,7)-x(4,8)).*(sin(x(4,9)+x(3,5))/sin(x(4,9)))).*x(4,4)./x(3,3))))./(x(4,1)*x(3,4)*cos(x(3,5))).*cos(x(3,8))+b.*(x(4,1)*x(3,4)*cos(x(3,5))).*(exp(-pi.*(x(2,:).*x(4,3)./x(3,2)).^2)-exp(-pi.*((x(2,:)-(x(3,7)-x(4,7)-x(4,8))).*x(4,3)./x(3,3)).^2)).*(x(4,1)*x(3,4)*cos(x(3,5))/2.*(erf(sqrt(pi).*x(1,:).*x(4,2)./x(3,1))-erf(sqrt(pi).*(x(1,:)-(x(3,6)-x(4,5)-x(4,6))).*x(4,2)./x(3,1))))./(x(4,1)*x(3,4)*cos(x(3,5))).*sin(x(3,8))','b','x');
[b,norm,res,ef,out,lam,jac]=lsqcurvefit(f,b,[x;y;z;w],u,[],[],options);
特别注意的是在这两个函数中都调用了LevenbergMarquardt(L-M)极小化方法，以避免出现局部解和病态问题。
4.3 预计参数的实例验证
本文以皖北某矿井1013工作面为例，结合地质采矿条件和观测点实测数据，利用MATLAB提供的函数求取地表移动与变形预计参数。此工作面走向长度575m，倾向长度150m，下山采深为415m，上山采深为361m，平均采深为388m，煤层平均倾角10o，煤层平均厚度3.1m。在该工作面上方地表布置两条观测线(走向和斜向)，由于受实地地形的限制，走向设置大半条观测线，斜向观测线观测点布设在水渠上，走向观测线与斜向观测线的夹角为44o。地表移动观测站布设如图2所示,工作面实测数据整理如图3所示。
预计的走向观测线下沉值和水平移动值的中误差分别为60.01mm和21.4mm，分别占最大值的2.6%和2.5%，预计的斜向观测线下沉值和水平移动值的中误差分别为73.3mm和34.4mm，分别占最大值的3.1%和4.5%。
5 结论
通过预计的拟合曲线与实测数据的可视化比较，可知只要充分利用MATLAB工具箱函数，结合任意工作面监测点实测数据及实际的地质采矿条件，就能够拟合出较为准确的地表移动变形预计参数，极为方便，但在此需特别指出，拟合时初始化的参数非常重要，因此要充分利用周围矿区的资料。
参考文献：
[1] 邹友峰,邓喀中,马伟民.矿山开采沉陷工程[M].徐州:中国矿业大学出版社,2003
[2] 吴侃,周鸣.矿山开采沉陷预测预报系统[M].徐州:中国矿业大学出版社,2000
[3] 康建荣,王金庄,温泽民.任意形多工作面多线段开采沉陷预计系统(MSPS)[J]矿山测量,2000,1(3):24~27
[4] 朱刘娟,陈俊杰,邹友峰.任意形状工作面开采地表移动变形预计的算法实现[J].辽宁工程技术大学学报,2005,24(3):337~339
[5] 陈良松,汪青松,杨永国. 巨厚松散层条件下地表移动变形特征.煤矿安全,2011 42(8)
[6] 王明柱,郭广礼,王磊,王彬.基于岭估计的概率积分法预计参数的求取.煤矿开采, 2012 17(2)
[7] 周万茂,张华兴,仲惟林.多工作面开采影响下求取地表移动参数方法研究.煤炭学报.2000 25(3)
[8] 梁明,王成绪.厚黄土覆盖山区开采沉陷预计.煤田地质与勘探.2001 29(2)
[9] 朱广轶,王玮,刘晓群.开采与地表沉陷分析的可视化软件系统.煤炭工程.2006 (8)