http://www.paper.edu.cn
基于MATLAB的开采沉陷预计系统1
李培现，谭志祥，齐公玉，顾伟，王磊
中国矿业大学环境与测绘学院，江苏徐州（221008） 江苏省资源环境信息工程重点实验室，江苏徐州（221008）
E-mail：pxlicumt@126.com
摘 要：本文介绍了利用MATLAB开发开采沉陷预计系统的实现方法。主要介绍了开采沉陷预计的数学模型，系统结构、功能、MATLAB在开采沉陷预计分析中的应用及优势。最后通过实例验证了程序的正确性，预计结果和手算结果吻合，求参数部分观测值和拟合值基本一致，证明了该系统计算结果是可靠的。系统利用MATLAB的图形工具实现了开采沉陷预计的可视化，与SURFER、GIS等有良好接口，便于后处理和分析。 关键词：开采沉陷；MATLAB；概率积分法；预计参数 中图分类号：TD73
据国有重点煤矿的不完全统计资料，目前我国“三下”压煤约137.9亿吨[1]，部分矿区“三下”压煤量十分巨大，已经严重制约着矿区的发展，如徐州矿务集团2006年底“三下”压煤量约3.1亿吨，占总可采储量的62%左右。如何正确地对一个计划进行的开采进行岩层和地表沉陷预计和有效表达对指导“三下”开采实践有重要的作用。目前，常用的开采沉陷预计软件大多使用高级编程语言VB、VC++编写[2~4]，编程效率较低，且大多程序在预计结果的可视化表达和处理方面存在很大的不足。
MATLAB是美国MathWorks公司研制开发的软件产品，它是当今世界上使用最为广泛的科学计算软件之一[5~7]。MATLAB具有强大、丰富的内置函数和图形显示功能，另外MATLAB还提供了灵活而又强大的接口技术和图形用户界面。因此MATLAB将数据结构、编程特性以及图形用户完美的结合到一起。本文研究了利用MATLAB的强大功能实现矿山开采沉陷的预计、预计结果的可视化表达以及地表移动实测数据处理的方法。
1. 开采沉陷预计的数学模型
概率积分法是基于随机介质理论的一种预计方法。因其算法简单、结果可靠，目前已成为我国较成熟的、应用最为广泛的预计方法之一[8~10]。因此，采用概率积分法编制了开采沉陷预计系统。概率积分法数学模型为：
① 地表任意点（x，y）的下沉值W（x，y）
W(x,y)=
其中
1
W0(x)W0(y)W0
（1）
W0(x)=W0(y)=
W0ππ[erf(x)−erf((x−l)] （2） 2rrW0ππ[erf(y)−erf((y−L)]2r1r2
（3）
（4）
W0=mqcosα 1本课题得到国家自然科学基金资助项目（40772191）；“十一五”国家科技支撑计划重点项目（2006BAC09B01）的资助。
- 1 -
http://www.paper.edu.cn
式中W0—充分采动条件下地表最大下沉值；m—采厚；q—下沉系数；α—煤层倾角；W0（x）、W0（y）—走向和倾向有限开采时主断面地表下沉值；l、L—走向和倾向有限开采时的计算长度（考虑拐点偏距后的长度）；r、r1、r2—分别为走向、下山、上山的主要
erf(x)=
影响半径；erf（x）—高斯误差函数，
2
e∫π0
x
−x2
dx
。
② 地表任意点（x，y）沿ϕ方向的倾斜值、曲率值、水平移动和水平变形值分别为：
i(x,y,ϕ)=
K(x,y,ϕ)=
10
[i(x)W0(y)cosϕ+i0(y)W0(x)sinϕ]W0
（5）
1
[K0(x)W0(y)cos2ϕ+W0(x)K0(y)sin2ϕ+i0(x)i0(y)sin2ϕW0
（6）
U(x,y,ϕ)=
1
[U0(x)W0(y)cosϕ+U0(y)W0(x)sinϕ]W0
（7）
ε(x,y,ϕ)=
102{ε(x)W0(y)cosϕ+ε0(y)W0(x)sin2ϕ+[U0(x)i0(y)+U0(y)i0(x)]sinϕcosϕ}W0
（8）
0
0
式中：i(x,y,ϕ)，i(x)，i(y)—（x，y）点沿ϕ方向、充分采动时（x，y）点对应
00
K(x,y,ϕ)K(x)K(y)—（x，y）点沿ϕ方向、走向断面及倾向主断面的倾斜值；，，00
U(x,y,ϕ)U(x)U充分采动时（x，y）点对应走向主断面及倾向主断面的曲率值；，，(y)—
（x，y）点沿ϕ方向、充分采动时（x，y）点对应走向断面及倾向主断面的水平移动值；
ε(x,y,ϕ)，ε0(x)，ε0(y)—（x，y）点沿ϕ方向、充分采动时（x，y）点对应走向断面
及倾向主断面的水平变形值
③ 多个工作面开采引起的下沉和变形预计考虑所有工作面对地表点影响的总和，采用叠加的方法处理。
2. 系统结构及功能
按照结构化和模块化的设计思想，将开采沉陷系统分为以下6个功能模块（如图1）。各模块均编写了相应函数和M文件。
开采沉陷预计系统预计数变形预计图形编据读入辑输出预计结参数求取果输出
图1 系统结构图 Figure 1 Structure of the System
各个模块功能如下：
（1） 预计数据读入模块：实现预计概率积分法参数和工作面的直接读入。预计计算前
- 2 -
http://www.paper.edu.cn
需提前准备预计的数据文件。
（2） 变形预计模块：利用概率积分法预计以下内容上的下沉和变形值： ①、矩形区域范围内； ②、主断面上； ③、任意离散点；
（3） 图形编辑输出模块：实现预计下沉及变形等值线二维图形的自动绘制及输出，下沉盆地三维效果图的自动绘制和输出。
（4） 预计结果输出：输出Surfer、GIS等的接口数据文件，以便与其他相关软件结合进行工程计算和分析。
（5） 参数求取：根据观测站实测资料，利用模矢法[2]求取概率积分法预计参数。
3. MATLAB在开采沉陷预计系统中的关键技术
（1） 高斯误差函数的计算：概率积分法中用到了高斯误差函数， MATLAB中包含了对高斯误差函数进行计算的方法和函数，命令erf（x）等价于计算高斯误差函数的值，即：
erf(x)=
2
π∫
x
0
e−xdx
（9）
2
（2） 网格的自动剖分：在一定区域的预计中，往往要根据一定的规则将预计区域划分为大量的网格点以在后续处理中进行等值线或者三维曲面的绘制。MATLAB中meshgrid命令的作用是将给定的区域按照一定的方式划分为平面网格，该网格可以用来绘制二维的等值线或者三维的曲面图。对于预计区域的网格剖分只需要使用meshgrid命令即可以实现。
（3） 绘图功能：MATLAB的数据分析和可视化功能可以表达出数据的二维、三维和四维的图形，这些为开采沉陷预计的可视化表达提供了很大的方便。
（4） 坐标转换：程序中坐标转换采用三个已知值：工作面左下山角点坐标（X0，Y0）和走向方位角α值来转换到概率积分法的计算坐标系。令大地坐标系（左手坐标系）或者其他坐标系（左手坐标系）中坐标为（X，Y）转换到计算坐标系（右手坐标系）中为（Xj，
Yj）。转换公式为：
⎡Xj⎤⎡cosα⎢⎥=⎢⎢⎣Yj⎥⎦⎣sinαsinα⎤
[X−X0
−cosα⎥⎦
Y−Y0] （10）
（5） 模矢法求参中误差函数的求取：程序中利用模矢法求取概率积分法参数提供了三种误差函数分别为：根据下沉观测值求参，根据水平移动值求参和根据下沉和水平移动联合求参数,即分别为：
2
[VV]=∑[(W(B)−W(实）]
（11） 2
[VV]=∑[(U(B)−U(实）]
（12）
22
[VV]=∑[(W(B)−W(实）]+∑[(U(B)−U(实）]
（13）
式中 [VV]—误差函数
W（B）、W（实）—概率积分法计算下沉值和实测下沉值 U（B）、U（实）—概率积分法计算水平移动值和实测水平移动值
误差函数分别为误差向量2—范数的平方，即（11）、（12）、（13）三式可等价的表示为：
- 3 -
http://www.paper.edu.cn
[VV]=
W(B)−W(实）2U(B)−U(实）2
2
2
2
（14） （15）
2
[VV]=
[VV]=
W(B)−W(实）2+U（B）—U（实）2
（16）
MATLAB含有丰富的矩阵运算函数，使用norm命令可以求取向量2—范数，这也使求参过程中误差函数的计算变得简单而方便。
以上所列开采沉陷预计和参数求取过程中的关键技术，在传统的高级编程语言中需要编写和调试大量程序才能完成，而在MATLAB中只需几个简单的命令。这是使用MATLAB编程简单、高效的体现，也是MATLAB在开采沉陷预计程序编制中的优势。
4. 程序正确性验证
4.1 预计部分程序正确性验证
采用参考文献[10]187页第四题中参数进行求解，概率积分参数为：下沉系数0.7，主要拐点偏距10m，水平移动系数0.3，开采影响传播角90°；工作面走向长150m，影响角正切2.0，
；工作面走倾向长100m，采深100m，采厚1m，水平煤层；工作面左下角坐标为（0，0）
向方位角为90°；预计点坐标为（100，0），预计方向为135°。表1给出了手算和程序计算结果，计算最大的相对误差为0.22%。图2～图5给出了由MATLAB直接生成的地表下沉、倾斜、曲率等值线图和下沉盆地三维效果图。由以上计算结果可见，开发的地表移动变形预计程序是正确的。
表 1 预计结果对比表 Tab.1 compare for the result of predict
计算手算程序
下沉倾斜曲率水平移水平变形
66.58066.434
5.385.37-0.00
0.31350.3133-0.000
80.77524.7038
80.686 4.7002 -0.02 -0.0031
差值 -0.145
图2 下沉等值线图 图3 倾斜等值线图 Fig.2 Subsidence contour Fig.3 Incline contour
- 4 -
http://www.paper.edu.cn
图4 曲率等值线图 图5 下沉三维图
Fig.4 curvature contour Fig.5 three-dimensional figure of subsidence
4.2 求取概率积分法参数部分程序正确性验证
根据参考文献[2]148页例子测试程序求参数的效果。工作面采厚m＝2.40米，煤层倾角
α＝11O，工作面倾向长D1＝180米，走向长D3＝250米，下山边界采深H1＝537米，上山边界采深H2＝484米，观测站实测下沉与水平移动值见表2[2]。
表2 实测下沉与水平移动值
Tab.2 observation of subsidence and horizontal displacement
点号 X坐标 Y坐标 下沉值 水平移动值 移动方向
1 -125 29 2 -106 35 3 -100 36 4 -67 48 5 -16 65 6 31 81
104 118 151 154 168 166 300 232 506 2 682 200
19 19 19 20 20 20 22 24 28 30 31 31
7 84 101 737 29 8 112 112 686 -69 10 167 137 11 194 150 12 220 165 13 228 168
456 -185 324 -184 207 -153 180 -142
9 140 124 585 -145 25
14 239 174 143 -123 31
利用本系统和表2中的观测数据求参数，结果和文献[2]中对比情况见表3。
表3 文献结果与程序计算结果对比表 Tab.3 Result compare of literature and programe
参数 q b tanβ S1 S2 S3 θ 文献值 0.51 0.25 2.49 -50 10 24 85 求取值 0.54 0.25 2.40 -38 10 28 85
- 5 -
http://www.paper.edu.cn
利用求得的参数对观测点进行预计。下沉曲线拟合的效果见图6。水平移动曲线拟合效果见图7。由图可知，利用本预计系统，根据观测站实测数据求取开采沉陷概率积分法参数的结果是正确可靠的。
参数拟合下沉值实测下沉值123456710111213143002001000-100-200123参数拟合水平移动值实测水平移动值020040060080045671011121314 -300图6 下沉值拟合图 图7 水平移动值拟合图
Fig.6 data fit of subsidence Fig.7 data fit of horizontal displacement
5. 结论
（1）基于概率积分法基本原理，利用MATLAB开发了开采沉陷预计系统，理论基础可靠，可以满足工程的需要。
（2）充分发挥了MATLAB的在绘图和数值计算方面的优势，生成地表移动变形的曲线图、等值线图、下沉盆地三维视图。大大减少了软件开发的工作量。
（3）可以对单一工作面的观测站数据进行处理，求取地表移动变形的概率积分法预计参数。
（4）输出了SURFER、GIS的接口数据文件，易于与其他软件结合进行后处理和分析。 （5）坐标转换由程序自动完成易于与大地坐标系或者坐标系下的矿图结合。
参考文献
[1] 谭志祥,邓喀中.建筑物下采煤理论与实践[M]. 徐州:中国矿业大学出版社,2007 [2] 吴侃,周鸣. 矿区沉陷预测预报系统[M]. 徐州:中国矿业大学出版社,1999.
[3] 朱庆伟,郭达志. 煤矿区地表沉陷及其可视化新方法[J]. 矿业安全与环保,2006（6）:10-12
[4] 赵林林,佟光成,申灵玲. 基于ArcObjects的地表沉陷信息可视化技术研究与实现[J].中国矿业,2007（11）:33-35
[5] 董维国编著.深入浅出MATLAB7.X混合编程[M].北京:机械工业出版社,2005.10 [6] 孙祥,徐流美,吴清编著.MATLAB7.0基础教程[M].北京:清华大学出版社,2005.5 [7] 求实科技.MATLAB7.0从入门到精通[M].北京:人民邮电出版社,2006.3
[8] 国家煤炭工业局. 建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程[M]. 北京:煤炭工业出版社. 2000.
[9] 邹友峰,邓喀中,马伟民. 矿山开采沉陷工程[M]. 徐州:中国矿业大学出版社.2003. [10] 何国清,杨伦等. 矿山开采沉陷学[M]. 徐州:中国矿业出版社.1991
- 6 -
http://www.paper.edu.cn
Mining Subsidence Prediction System Based on MATLAB
Li Pei-xian, Tan Zhi-xiang, Qi Gong-yu, Gu Wei, Wang Lei
School of Environmental and Spatial Information in China University of Mining & Technology,
Jiangsu Xuzhou (221008)
Key Laboratory of Resources and Environment Information Engineering in Jiangsu Province,
Jiangsu Xuzhou (221008) Abstract
This paper introduced the method to program Mining Subsidence Predict and Analysis System based on MATLAB. Mainly on theory of Mining Subsidence, structure and function of the system , how Matlab used and advantage in Mining Subsidence Predict and Analysis. An example was used to validate the programe at last, the results of predict tallies the result computed with hand, the result of observation and computed with new parameters were consistent in part of parameters computed. Those test proved the results of the system are credibility. The system realized visualization of Mining Subsidence Predict, have good interface with SUEFER and GIS, covenient for after process and analysis.
Keyword：Mining Subsidence, MATLAB, Probability integration Method, Parameters
作者简介：李培现(1983-)，男，山东巨野人，硕士研究生，从事矿山开采沉陷与数字矿山方面的研究工作。
- 7 -
