基于ArcGIS Engine的开采沉陷数据处理和可视化

李青青，徐良骥，翟应峰

（1.安徽理工大学 测绘学院，安徽 淮南232001）

基于ArcGIS Engine的开采沉陷数据处理和可视化

李青青1，徐良骥1，翟应峰1

（1.安徽理工大学 测绘学院，安徽 淮南232001）

调用工作面内的参数，利用ArcGIS Engine通过概率积分法计算得出工作面内任意点的下沉值，并根据得出的沉陷值生成图像。这样就可以比较直观地将地表的变形情况展现出来，既能输出具体变形数值，又能够以图像的形式表达。

Arc Objects;沉陷;参数;概率积分法;图像

GIS的发展与普及为矿山开采沉陷控制提供了一种数据处理与可视化的方法与途径[1]。对于采区的变形监测数据进行处理，根据主要影响参数计算出沉陷值，包括数据处理和可视化2个主要部分。此文基于ArcObjects进行开采沉陷数据处理与可视化的表达。因为地表下沉比较抽象，地下开采后对地表的影响需要三维图像对其进行逼真的描绘和定量化的表达，利用已知的数据对其进行三维可视化表达，给地下开采后地表沉陷比较直观形象的图像[2,3]，便于直观的观察和研究。

1 关键技术简介

1.1 ArcGIS Engine技术

ArcGIS Engine（AE）是用于构建定制应用的一个完整的嵌入式GIS组件库，它可以为用户提供针对GIS解决方案的定制应用[4]。利用 ArcGIS Engine ，开发者可以将 ArcGIS 功能集成到一些应用软件中，但它并不是一个最终产品，而是一组ArcObjects 软件开发包[5]。

ArcGIS Engine可以在Windows、UNIX和Linux桌面上运行并支持一系列的应用软件开发环境，包括Visual Basic 6，Microsoft Visual Studio NET和各种Java开发环境如表1[6]。

表1 ArcGIS Engine支持一系列操作系统平台和编程语言

此文主要是利用了ArcGIS Engine提供的几种控件——MapControl、ToolbarControl、PageLayOutControl，以及用于三维场景显示的SceneControl等进行开发，并且对于开发出的应用程序赋予了一定的license，可以单独运行于没有任何ArcGIS程序的微机上。

1.2 Visual C#.NET简介

Visual C#．NET是一套综合工具集，用于为Microsoft Windows和Web创建XML Web服务和基于Microsoft．NET的应用程序。Visual C#．NET是基于强大的C++传统语言创建的。C#是一种现代化的面向对象的程序开发语言，它使程序员能够在微软．NET平台上快捷开发种类丰富的应用程序。．NET平台提供了大量的工具盒服务，能够最大限度地挖掘和使用计算及通信能力[7]。

系统以微软的Visual Studio．NET为开发环境，它是一个功能强大、高效并且可扩展的编程环境[8]。在语言上主要基于C#编程为主体程序并且调用其他功能模块实现系统的功能。

2 概率积分法简介

概率积分法是因其所用的移动和变形预计公式中含有概率积分（或其导数）而得名。由于这种方法的基础是随机介质理论，所以又叫随机介质理论[9]。在二维情况下，若开采单元的横坐标为s,地表任意点A的坐标为x,则此单元开采引起的A点的下沉为：

若煤层是水平的，煤层坐标系（t，s）和地面坐标（x，y）的水平投影重合，则单元B（坐标为（s，t））的开采引起地表任意点a（坐标为（x，y））的下沉W0（x，y）为[10]：

3 系统设计思路

程序是根据概率积分法思想来实现的。从统计观点看，概率积分法认为任意开采条件下都可以把整个开采分解为许多或无线多个微小单元的开采，整个开采对地表的影响等于所有单元影响的总和，而可以从单元开采入手，研究下沉盆地的方程式。主程序流程如图1所示。

图1 主程序流程图

4 系统功能与实现

此文基于AE编程，调用工作面参数，根据概率积分法原理，计算得出工作面内任意点的下沉值，并根据沉陷值生成图像[11]。计算中主要用到MapControl、TOCControl、ToolbarControl、SceneControl4个控件。需要说明的是，每个AE程序都要加一个LicenseControl控件。

操作基本上按图1主程序流程图进行，主要分为以下几个步骤：

第1步，数据处理。点击“生成下沉值”按钮，输入工作面参数，确定后输入保存路径，保存处理后的数据。

第2步，数据格式转换，将txt文本转换为Shapefile格式的文件。

第3步，加载数据，生成点状图像。

第4步，生成TIN，并在MapControl与SceneControl中同步显示。

4.1 数据处理

首先，需要收集采区的部分参数。本设计用到的采区参数主要包括：平均开采厚度、下沉系数、煤层倾角弧度、走向岩层厚度、上山覆岩厚度、下山覆岩厚度、走向长度、倾向长度、左边界拐点偏移距、右边界拐点偏移距、上边界拐点偏移距、下边界拐点偏移距、最大影响角正切值，这些数值收集和整理后可以在系统运行的时候输入到沉陷预计参数界面里面。

最大影响半径、工作面的长度和宽度、最大下沉值都可以根据以上这些参数通过概率积分法计算得出，然后计算工作面内任意点的下沉值并输出最终的下沉数据。下面给出计算任意点下沉值的方法。

在开采水平煤层时，设煤层采出厚度为m，开采面积为F，如图2所示，则开采单元积dF引起地表点的下沉将为DW0=q·mf（S）dF，当面积F全部采出后，地表点A的总下沉值为WA[12]。

式中，q为下沉系数，主要决定于顶板管理方法；f（S）为单位面积dF至A点距离S的函数，这个函数一般称为影响函数。

将f（S）、dF分成沿x、y 两个方向来考虑：

图 2 水平煤层

根据《几何理论》研究地表变形的规律，认为地下开采对地表点下沉的影响服从于正态分布，因此影响曲线函数为一高斯曲线，即

所以，

系统运行后，第一步操作是进行数据处理，点击“生成下沉值”，输入工作面参数后，点击“确定”按钮，在弹出的对话框中输入保存路径，点击“确定”后，系统会自动弹出“完成”对话框，点击并完成第一步操作。数据处理后将以txt文本格式输出。

4.2 数据格式转换

第一步处理完数据生成下沉值后，以txt文本格式输出。此时，需要转换数据格式才能进行下一步操作，类似于ArcGIS中数据格式之间的转换。设计中主要是将txt文本格式转换成Shapefile格式，然后再将其加载进程序中。

从txt到Shapefile转换主要有2种方式。其一主要是在ArcView中进行。具体步骤是file--extentions--Generate to shape，可以加载相应的模块，然后点击view--add Generate Theme对话框中加载，完成转换。此文采用另一种方式，利用ArcToolbox/Samples/Data Management/Features中 的2个工具Create Features From Text File（导入）和Write Features To Text File（导出）完成txt文件格式的要素导入导出。程序流程如图3所示。

图3 txt转Shapefile流程图

格式转换的部分代码如下：

OpenFileDialog pOpenFileDialog = new OpenFileDialog();

pOpenFileDialog.Filter = "文本文件|*.txt";

pOpenFileDialog.InitialDirectory = Application. ExecutablePath;

if (pOpenFileDialog.ShowDialog() != DialogResult. OK) return;

m_txtFilePath = pOpenFileDialog.FileName.ToString();

SaveFileDialog pSaveFileDialog = new SaveFileDialog();

pSaveFileDialog.Filter = "shp文件|*.shp";

pSaveFileDialog.InitialDirectory = Application. ExecutablePath;

if (pSaveFileDialog.ShowDialog() == DialogResult. OK)

……

系统运行时，点击“txt转换为shp”按钮，在弹出的对话框中选择第一步生成下沉值后输出的txt文本，选中并确定后，选择保存路径将转换格式后的数据命名、保存起来。格式转换完成以后，需要检查其最终的数据是否正确，如果发现空文本或空数据集，说明转换不正确。还可以通过加载格式转换后的数据文件到系统中，进一步验证。

4.3 图像显示

数据转换为Shapefile格式后，就可以将其加载进来。

添加shp文件的代码如下:

//添加shp文件到Map

private void AddData_Click(object sender, EventArgs e)

{

arclss.Addshp();

this.axMapControl1.Map.AddLayer(arclss.PLayer as ILayer);

this.axMapControl1.Refresh();

点击“加载数据”，选中格式转换后的文件，点击“确定”，成功添加shp文件后，图层树控件——TOCControl中会有一个shp文件的图层，它可以显示出文件的要素为点状要素，而右边地图控件——MapControl中显示的图像则为本次设计选取的所有采样点。

加载数据成功之后，就可以在工具栏控件——ToolbarControl上面进行操作，如中心放大、中心缩小、漫游、全图、选中位置缩小和放大等。

通过这些点状图像，并不能直观地获取其他信息，因此，接下来需要通过其生成TIN，然后加载TIN进来，在MapControl与SceneControl中分别显示二维和三维的图像。从散点到连续表面，需要引进数字高程模型的概念。

4.4 数字高程模型

数字高程模型是通过有限的地形高程数据实现对地形曲面的数字化模拟。最为常见的表示方法主要是规则格网——GRID和不规则三角网——TIN，如图4所示。规则格网的数据组织类似于图像栅格数据，而它的每个像元的值是高程值，即格网DEM是一种高程矩阵。不规则三角网直接利用原始的采样点进行地表重建，由连续的相互连接的三角面组成，三角面的形状和大小取决于不规则分布的观测点的密度和位置。

图4 GRID和TIN

本设计主要通过TIN来建立地形模型。在格式转换后，运行系统，点击“三维显示”按钮，在弹出的对话框中选中转换格式后的shp文件并确定，生成并加载TIN,如图5所示。

图5 生成TIN图像

生成TIN的部分代码如下：

pOpenFileDialog.Filter = "shp文件|*.shp";

if (pOpenFileDialog.ShowDialog() == DialogResult. OK)

{

string tempPath = pOpenFileDialog.FileName. ToString();

stringtemp1=System.IO.Path.GetDirectoryName (tempPath);

string tempname =System.IO.Path.GetFileName (tempPath);

IWorkspaceFactoryp WorkspaceFact = newShapefile WorkspaceFactoryClass();

IWorkspacepWorkspace=pWorkspaceFact. OpenFromFile(temp1, 0);

IFeatureWorkspacepFeatureWorkspace=(IFeatureWor kspace)pWorkspace; IFeatureClasspFeatClass=pFeature Workspace.OpenFeatureClass(tempname);

IFeatureLayer pFLayer = new FeatureLayer();

pFLayer.FeatureClass = pFeatClass;

IFields pFields = pFeatClass.Fields;

IFieldpHeighField=pFields.get_Field(pFields. FindField("S"));

IGeoDataset pGDS = pFeatClass as IGeoDataset;

IEnvelope pEnv = pGDS.Extent;

pEnv.SpatialReference = pGDS.SpatialReference;

ITinEdit pTinEdit = new TinClass();

pTinEdit.InitNew(pEnv);

object obj = Type.Missing;

pTinEdit.AddFromFeatureClass(pFeatClass, null, pHeighField, null, esriTinSurfaceType.esriTinMassPoint, ref obj);

pTinEdit.SaveAs("D:\daaa", ref obj);

pTinEdit.StopEditing(false);

在生成TIN并加载进来后，于MapControl与SceneControl中分别显示二维和三维图像。需要考虑二者的同步显示和更新问题，即同步缩小和放大等。以中心放大为例，部分代码如下：

IEnvelope penv = this.axMapControl1.Extent; penv.Expand(0.8, 0.8, true);

this.axMapControl1.Extent = penv;

IActiveView pActiveView = this.axMapControl1.Map as IActiveView;//得到接口

IEnvelope pEnvelope = (IEnvelope)pActiveView. Extent;//获得显示范围

pEnvelope.Expand(0.8, 0.8, true);

pActiveView.Extent = pEnvelope;

pActiveView.Refresh(); //刷新

……

ISceneViewer pSceneViewer = this.axSceneControl1. SceneGraph.ActiveViewer as ISceneViewer;

pSceneViewer.Redraw(false);//重绘现场

4.5 三维场景

生成TIN图像后，在MapControl中的二维显示如图5，还可以进一步通过分层设色来使图像中地貌高程分布及其相互对照鲜明，如图6所示。但是，通过2 个图我们仍然不能够很直观地看出地表的变形，需要通过场景控件——SceneControl来显示三维场景，比较直观地表达地表的变形情况，显示出来的图像也是整个设计最终的可视化结果，如图7所示。

图6 分层设色

图7 地表变形三维场景

我们可以通过ToolbarControl上的工具来调整图像，如点击“漫游”，按钮在图中可以观察任一方向的地表变形。

5 结 语

1）本文以概率积分法为基础，利用C#语言开发出开采沉陷预计程序核心程序，通过调用工作面参数预计地表下沉盆地内任一点的下沉值，并输出数据。

2）在设计中完成了数据格式的转换，将输出数据的txt文本格式转换为Shapefile格式，并利用AE将转换格式后的数据以图形输出，这样就可以比较直观地观察工作面内的点分布。

3）完成了TIN的构建，从而实现了开采沉陷可视化，可以直观地观察到工作面内的地表变形情况。

4）实现了MapControl与SceneControl的同步显示和更新，并初步实现分层设色，使图像中地貌高程分布及其相互对照更加鲜明。

[1] 赵林林,佟光成,申灵玲,等． 基于ArcObjects的地表沉陷信息可视化技术研究与实现[J]．中国矿业,2007(11):33-35

[2] 王飞．基于MapGIS的开采沉陷预计分析系统[D]．焦作：河南理工大学，2010

[3] 路兵．基于C#的矿山开采沉陷预计可视化系统[D]．青岛：山东科技大学，2011

[4] 李青元，林宗坚，李成名．真三维 GIS 技术研究的现状与发展[J]．测绘科学，2000,25(2)：47-51

[5] 周新耿，刘芬．二维三维空间信息系统的结合和应用[J]．计算机与数字工程，2005,33(3):101-103

[6] 刘仁义，刘南．ArcGIS开发宝典——从入门到精通[M]．北京：科学出版社，2006

[7] 刘贵明,张卫国,兰小机． 基于ArcGIS Engine的矿山3D-GIS的设计与实现[J]． 金属矿山,2009(5):92-95

[8] 王玲，吴侃，奚新丽． 基于ArcGIS Engine的开采沉陷预计系统[J]． 测绘科学,2008,33(4):153-154

[9] 吴侃． 开采沉陷预计一体化方法[M]． 徐州：中国矿业大学出版社，1998

[10] 张应昌,孙亚峰,何尧． 基于概率积分法的开采沉陷预计分析[J]．山西建筑,2009(20):119-120

[11] 白力改，谭志祥，柳聪亮，等．基于ArcGIS的开采沉陷预计分析可视化方法[J]．金属矿山,2010(5):140-142

[12] 翟厥成,洪加明．地表移动盆地内任意点下沉值的预计[J]．矿山测量,1978(2):35-43

[13] Demers M N． 地理信息系统基本原理 [M]．北京:电子工业出版社，1999

[14] 张新宇,陈殿友,杨天行． Arc/View中Shapefile数据结构剖析及格式转换[J]． 吉林大学学报:信息科学版,2003(4):399-402

P208

B

1672-4623（2015）02-0126-04

10.3969/j.issn.1672-4623.2015.02.045

李青青，硕士，研究方向为地图制图与地理信息工程。

2014-03-12。

项目来源：安徽省自然科学基金资助项目（1208085QE91）；安徽省国土资源科技资助项目（2012-k-24）。