1、基于 MATLAB 的沉陷预计参数求取摘要：本文利用 MATLAB 工具箱函数，设计了求取地表移动变形预计参数的程序，结合工作面任意点的观测数据和地质采矿条件，能够很好的求取地表移动变形预计参数。本文通过实例，对所拟合的参数进行检验，证明了此方法的可行性。 关键字：MATLAB；非线性最小二乘法；概率积分法；预计参数 Abstract: This paper designs the program to calculate the surface movement parameters by using MATLAB toolbox functions. By combining with t
2、he arbitrary point observation data and the conditions of geological mining, The program can obtain the surface prediction parameters. In this paper, the estimated parameters are inspected, the feasibility of this method is proved by the results of experiments. Keywords: MATLAB; nonlinear least squa
3、re method; probability integral method; predicting parameters 中图分类号：TU441+.7 文献标识码：A 文章编号：2095-2104（2013） 1 引言 MATLAB 有美国 MathWorks 公司推出，其与 Mathematica、Maple 并称为三大数学软件，具有优秀的数值计算能力和卓越的数据可视化功能。 MATLAB 自问世以来，已经发展称为适合多学科、多种工作平台的功能强大的大型软件，广泛应用于大学科研、工程计算等领域，尤其在工程界。概率积分法通过积分而得到下沉曲线表达式，因其所用的移动和变形预计公式中含有概率积分
4、而得名。这种方法是将矿山岩层移动作为一种服从统计规律的随机现象来讨论。因此，此种方法是以随即介质理论为基础的一种预计方法。概率积分法自 20 世纪 60 年代引入我国以来，在许多矿山得到广泛使用，并不断发展，目前是我国进行地表移动变形预计主要的方法之一。地表移动变形预计参数对开采沉陷预计的准确性具有重要意义，因而本文以概率积分法数学模型为基础，借助 MATLAB 平台，编写求取地表移动变形参数的程序。 2 概率积分法求参的数学模型 利用概率积分法关于任意点的地表下沉预计公式作为求取预计参数的基本数学模型，即： 地表移动观测站任何一个测点的下沉实测值 W 全能表达成自变量(x,y)(观测点的水平
5、坐标)和概率积分法预计参数q、tg、s3、s4、tg1、tg2、s1、s2、 的函数，见下式： 在求取预计参数 b 时，则用到任意点的地表水平移动公式： 假设公式(3)里的参数q、tg、s3、s4、tg1、tg2、s1、s2、ctg，已经求出，那么上面的数学公式中就只剩下水平移动系数 b 这一待求参数。 3 Matlab 曲线拟合求参函数 Matlab 强大的绘图和用户自定义函数等功能，在数据处理领域被广泛应用，其中当然包括对实验、实测数据的拟合工作。鉴于本文研究内容，着重介绍 Matlab 中非线性曲线拟合（非线性最小二乘拟合）函数。 Matlab 中提供了两个求非线性最小二乘拟合的函数：l
6、sqcurvefit和 lsqnonlin。lsqcurvefit 的定义如下： 已知拟合的数据点为：xdata 和 ydata，lsqcurvefit 用以求含参量variable(向量)的向量值函数 F(x,xdata)中的参变量 x(向量)，使满足下式： (4) 完整调用格式为： 式中：x 是拟合的参数最优解，norm 是最优解的误差的平方和，res为误差向量，ef 是程序结束时的状态指示：ef0 表示结果收敛，ef=0 表示函数调用次数或迭代次数达到最大值（该值在 options 中指定） ，ef0表示不收敛，out 则包含数据的一个结构变量，包含实际的迭代次数等，lam 是上下界所对
7、应的 Lagrange 乘子，jac 是结果（x 点）处的雅可比矩阵，x0 为初始参数（缺省时程序自动取 x0=0）fun 是给出的目标函数，xdata 和 ydata 为拟合的数据，v1,v2 是拟合参数的上下界，可以用两个代替，options 是包含算法控制参数的结构设定（或显示）控制参数的命令为 Optimset 函数。 lsqnonlin 函数与 lsqcurvefit 函数的使用类似。 其完整的调用格式为 Matlab 可以很快速的直接读写 txt 格式和 Excel 格式的数据文件，所以给开采沉陷数据处理带来了方便，解决了求参和预计过程中繁多的数据手工输入和输出的问题。在 Matl
8、ab 中读写 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 地表移动变形
9、参数求取与检验 4.1 地表移动变形预计参数求取的算法流程 算法具体实现流程如图 1 所示： 图 1 求取参数的算法流程 4.2 求取沉陷预计参数的程序设计 本文主要用到的两个自定义函数 budgets 和 horizonal_budgets，函数 budgets 用于求取沉陷预计参数q、tg、s3、s4、tg1、tg2、s1、s2、，部分源代码如下： options=optimset(Display,iter,MaxFunEvals,28,MaxIter,214,TolX,1e-3,TolFun,1e-4,LargeScale,off,LevenbergMarquardt,on,LineSe
10、archType,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
11、(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,28,MaxIter,214,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
12、,:).*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
13、(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,j
14、ac=lsqcurvefit(f,b,x;y;z;w,u,options); 特别注意的是在这两个函数中都调用了 LevenbergMarquardt（L-M）极小化方法，以避免出现局部解和病态问题。 4.3 预计参数的实例验证 本文以皖北某矿井 1013 工作面为例，结合地质采矿条件和观测点实测数据，利用 MATLAB 提供的函数求取地表移动与变形预计参数。此工作面走向长度 575m，倾向长度 150m，下山采深为 415m，上山采深为361m，平均采深为 388m，煤层平均倾角 10o，煤层平均厚度 3.1m。在该工作面上方地表布置两条观测线(走向和斜向)，由于受实地地形的限制，走向设置大
15、半条观测线，斜向观测线观测点布设在水渠上，走向观测线与斜向观测线的夹角为 44o。地表移动观测站布设如图 2 所示,工作面实测数据整理如图 3 所示。 图 2 工作面地表移动观测站平面布置图 图 3 工作面实测数据汇总 利用本文设计的程序将以上数据计算后结果如下表 1 所示。 表 1 初始参数及拟合参数 注：拐点移动距负为内移，正为外移。 将上述求得的采后预计参数，分别预计走向观测线和斜向观测线测点的下沉值和水平移动值，进而与实测数据对比，结果如图 4 所示。 (a)走向下沉预计值与实测值对比 (b)走向水平移动预计值与实测值对比 (c)斜向下沉预计值与实测值对比(d)斜向水平移动预计值与实测
16、值对比 图 4 预计值与实测值对比 注：走向设计线因为与 X 轴近似水平，所以走向拟合曲线在二维图形中显示。斜向设计线与 X 轴不垂直，所以斜向拟合曲线在三维图形中显示，并利用 MATLAB 视角函数 view,使曲线所在的平面正对观察者，便于观察。 预计的走向观测线下沉值和水平移动值的中误差分别为 60.01mm 和21.4mm，分别占最大值的 2.6%和 2.5%，预计的斜向观测线下沉值和水平移动值的中误差分别为 73.3mm 和 34.4mm，分别占最大值的 3.1%和4.5%。 5 结论 通过预计的拟合曲线与实测数据的可视化比较，可知只要充分利用MATLAB 工具箱函数，结合任意工作面
17、监测点实测数据及实际的地质采矿条件，就能够拟合出较为准确的地表移动变形预计参数，极为方便，但在此需特别指出，拟合时初始化的参数非常重要，因此要充分利用周围矿区的资料。 参考文献： 1 邹友峰,邓喀中,马伟民.矿山开采沉陷工程M.徐州:中国矿业大学出版社,2003 2 吴侃,周鸣.矿山开采沉陷预测预报系统M.徐州:中国矿业大学出版社,2000 3 康建荣,王金庄,温泽民.任意形多工作面多线段开采沉陷预计系统(MSPS)J矿山测量,2000,1(3):2427 4 朱刘娟,陈俊杰,邹友峰.任意形状工作面开采地表移动变形预计的算法实现J.辽宁工程技术大学学报,2005,24(3):337339 5 陈良松,汪青松,杨永国. 巨厚松散层条件下地表移动变形特征.煤矿安全,2011 42（8） 6 王明柱,郭广礼,王磊,王彬.基于岭估计的概率积分法预计参数的求取.煤矿开采, 2012 17(2) 7 周万茂,张华兴,仲惟林.多工作面开采影响下求取地表移动参数方法研究.煤炭学报.2000 25（3） 8 梁明,王成绪.厚黄土覆盖山区开采沉陷预计.煤田地质与勘探.2001 29(2) 9 朱广轶,王玮,刘晓群.开采与地表沉陷分析的可视化软件系统.煤炭工程.2006 (8)
