基于GIS的开采沉陷灾害预测系统

2013-12-07郑朝治谭荣建陈晓军

浙江农业科学 2013年11期

郑朝治，谭荣建，陈晓军

(1.昆明理工大学国土资源工程学院，云南昆明 650093;2.云南省文山州国土资源局土地整理中心，云南文山 663099)

矿体被采出之后，采空区上覆岩体的应力平衡状态受到破坏，引起地表变形、移动的现象称为开采沉陷［1］。开采沉陷是一种作用力集中、危害巨大的地质灾害。它给矿区的生产、生活造成巨大的破坏，如铁路架空、道路裂缝、耕地减少、房屋倒塌等。尤其是在人口密集、经济发达的城镇及其附近，采矿与保护环境、建筑物安全的矛盾越来越突出。据不完全统计资料，目前我国“三下”压煤约137.9亿t［2］，部分矿区“三下”压煤量十分巨大，已经严重制约了矿区的发展。

开采沉陷预测的基本目的是对开采沉陷及其对环境产生的影响做出合理的评价，从而避免开采对地表建筑物等造成严重破坏，有助于对沉陷情况的控制和矿区项目的决策。

20世纪90年代以来，随着计算机技术的进一步发展，以及地理信息系统 (GIS)产业在全世界普及和发展，GIS突出的空间分析与三维可视化功能逐渐应用于开采沉陷问题的研究。通过GIS软件平台的二次开发和定制，可以较容易地实现预测和可视化的程序化，提高预测的速度和精度。作者对以C#为语言工具，以ArcGIS作为开发主要技术平台的开采沉陷灾害预测系统做出设计与研究。该预测系统对沉陷数据进行空间分析和可视化，有利于对地表移动和破坏的更直观和更深刻的理解。

1 开采沉陷预测的数学模型

研究开采沉陷作为一门学科始于20世纪20年代，国内外众多学者为此做出了许多尝试和努力。提出的预测方法和沉陷模型大致可以分为经验方法、剖面函数法、影响函数法、理论模型——物理模型。

1.1 概率积分法

概率积分法属于影响函数法，因其移动和变形预测公式中含有概率积分而得名。它认为上覆岩土层是由裂隙所切割的大量极小碎块体所组成，地下开采引起的岩土与表层的移动规律和随机介质模型中碎粒的移动规律宏观上相似。由于概率积分法算法简单可靠，易于计算机实现，目前已成为我国较为成熟、应用最广的煤矿开采沉陷预测方法。

1.2 塌陷形态和剖面特点

开采刚开始时，随着采空面积增大，移动盆地的面积及地表最大下沉值也增大，盆地呈尖底的“碗状”，此时地表的采动影响称为非充分采动，如图1中过程1和2。随着开采工作面的推进，当最大下沉值达到最大值并不再随着采空面积增大而增大，但盆底仍呈尖底“碗状”，此时，地表的采动影响称为充分采动，如图1中过程3。工作面再继续推进，下沉最大值不再增大，此时盆地呈平底的“盘状”，此时的采动影响称为超充分采动［3］，如图1中过程4。半无限开采，是指在断面方向采空区为无限长，它也是超充分采动的一个理想状态。

沉陷剖面的一侧如图2所示，其中r是主要影响半径，锐角β为主要影响角，主要影响角正切tg β=H/r，H为采深。B是计算时认为的煤壁位置，D为拐点，O为实际煤壁位置。O到B的距离就是拐点偏距S0的值，当B在采空区那一侧时S0为正。

图1 开采过程中塌陷盆地形态的变化

图2 沉陷剖面示意图

1.3 地表移动变形计算

概率积分法中任意点 (x，y)的下沉值 W(x，y)可用下式计算:

式中:m为开采厚度;q为下沉系数;a为煤层倾角;W0(x)W0(y)为走向l和倾向l有限开采时主断面地表下沉值;erf(x)是高斯误差函数:r、r1、r2分别为走向、下山、上山的主要影响半径。

地表任意点 (x，y)沿方向φ的倾斜值、曲率值、水平移动值和水平变形值分别为:

式中:i(x，y，φ)，i0(x)，i0(y)是沿 φ方向、充分采动时点 (x，y)对应走向断面及倾向主断面的倾斜值;K(x，y，φ)，K0(x)，K0(y)是沿φ方向、充分采动时点 (x，y)对应走向主断面及倾向主断面的曲率值;U(x，y，φ)，U0(x)，U0(y)是沿 φ方向、充分采动时点(x，y)对应走向断面及倾向主断面的水平移动值;ε (x，y，φ)，ε0(x)，ε0(y)是沿 φ方向、充分采动时点 (x，y)对应走向断面及倾向主断面的水平变形值。

1.4 塌陷参数性质和获取方法

从概率积分法模型得出:走向方向达到充分采动的条件时，工作面走向长=l+S1+S2≥2r+S1+S2，即l≥2r;拐点处的下沉值为Wmax/2;拐点处的倾斜和水平移动值达到最大;距拐点左右0.4r处的曲率和水平移动值达到最大。

可以根据以上性质通过实测来获取拐点偏距S0和影响半径r的值，也可以通过它们来验证数学模型的正确性。另外，一些现有的经验数据也可以使用，例如:我国拐点偏距经验值为0.05～0.3 H［1］。由于实际测量误差，可能会造成多种途径得出的参数数据不完全一致，根据实际情况选定合适的参数值，才能保证系统预测结果的精度。

2 系统设计

2.1 ArcEngine组件

ArcEngine组件是ArcObjects组件跨平台应用的核心集合，它提供了多种开发的接口，可以适应NET，Java，VB和C++等多种开发环境。开发者可以使用这些组件来开发定制GIS和地图应用。目前，ArcEngine技术是开发组件式GIS所采用的主流技术，它逐渐代替了组件式 GIS的其他技术，如MO、AO等。使用ArcEngine开发出来的系统不需要安装桌面版ArcGIS，只需安装ArcEngine运行环

境即可，这样既节省了开发成本，又节约了系统运行所消耗的计算机系统资源［4］。

2.2 系统结构和功能

经过需求分析后，对软件进行初步的整体设计。然后逐一实现数据录入、地表沉陷数据的计算、沉陷数据的图形显示、数据转换接口等各功能部分。调试运行通过以后，将它们集成组织起来，形成系统总体的功能模块 (图3)。这样逐一进行，最终完成系统的总体开发。

图3 系统的设计方案

参数输入模块。预测开采沉陷，首先需要确定采区的基本情况，主要包括:开采厚度m，下沉系数q，煤层倾角α，主要影响角正切tgβ，采深H，拐点偏距S0等。

沉陷数据计算模块。将基本参数传入C#编写的数学模型，计算出所需的高程、距离或者角度数据，并通过ArcEngine组件开发的模块生成三维空间数据。同时也要支持生成过程数据，以方便对系统计算过程更直观的观察与正确性检验。

基于数据库应用的空间分析模块。通过ArcSDE获得地表地形数据，区县行政区划分数据，房屋，道路，河流等已有数据，并与沉陷空间数据进行叠加分析，得到沉陷区地表变化情况，建筑物受损情况，受损设施分布情况［5］。

成果的输出模块。系统的最终成果是多样性的，包括兴趣点单点沉陷查询、任一点的移动形变的变化图、受灾区域的面数据等。还可以加入简单的经济损失模型以计算出开采塌陷损失。

3 系统开发的关键问题

3.1 数值的使用方法

计算模块根据开采沉陷模型计算得到的是一系列的数据点阵，其中包括横纵坐标与点的下沉、曲率、水平移动和变形等值。这些数据点都是离散的，不能构成所需要的等值线结果。但是由于塌陷模型的特殊性，得到的沉陷面一定是关于地理坐标连续且平滑的。所以，可以通过插值的方法用这些离散点构建出很近似于真实情况的曲面。插值方法的实质是在每2个已知点间插入多个点，然后将这些点用直线段依次连接起来［6］。当加入了大量内插点时，用直线段连接后就可近似看作是一条光滑曲线。插值方法有许多种，其中3次样条插值法已经广泛应用于绘制等高线等解决实际问题的多种工程软件中，有较高的可靠性和较成熟的算法。所以本系统采用3次样条作为模拟曲面的算法。图4为沉陷数据模块并插值后得到的沉陷盆地效果图。

图4 矩形开采区的沉陷效果

3.2 地物的查询分析

模拟预测矿区的受灾影响，需要首先建立该区域关于受灾参数的等值线，例如竖直位移，水平位移的等值线，从而得到已达到灾害程度的面状要素。然后，通过进一步空间查询才能得到所影响的占有一定面积的地面单位，例如房屋，烟囱，道路等。基于空间位置的查询，是根据要素与要素之间的空间关系进行的。它用到 ArcEngine中的SpatialFilter组件类，其中SpatialRel用于获取或者设置查询时的空间关系。设置空间关系为eariSpatialRelIntersects(相交)和 eariSpatialRel(包含)时，用受灾面对其他地面设施面进行查询，就可以得到受灾的地物。除此之外还有另一种受灾情况需要考虑，那就是管网类的设施受损，比如公路、沟渠。此时必须得出受损的部分，仅仅查询计算并不能满足需求，必须对其进行空间分析，具体方法为:用上面的方法查询得到管线类面或线，然后利用 ItopologicalOperator接口的Intersection方法得到与受灾区域的交集部分。

3.3 数据的叠加分析

由于空间数据三维问题的复杂性，现有的软件大多数停留在生成沉陷等值线这一步。为了使系统得到的模型更加直观地表达，还需要建立塌陷前后3D模型的对比。使用 ArcEngine中 Scene加载DEM数据，使用 ApplyTin3Dproperties应用3D属性到遥感栅格图像［7］，这样就可以把影像覆盖在建立好的高程模型上。

系统先将输入的采空区数据进行计算，得到离散的点阵列数据，进行等高线绘制等操作。但是把竖直下陷数据和地表数据叠加时，必须重新自动获取离散点信息。这样在3次样条模拟的曲面上进行离散点阵列的取样，则会造成精度的损失。所以本系统使用最初计算模块采用的点阵列，经过坐标转换后直接与地表模型的对应点高程相加，然后更新地表DEM数据。这样解决数据叠加的问题就更加准确合理。

以上方法同样适用于多个工作面引起沉陷的叠加，多个工作面在计算模块中使用相同的离散点阵采样，先把离散点阵上变形值相加，再与地面空间数据叠加以获得相对较好的精度。