用ArcGIS的Python脚本计算采煤沉陷区土地损毁程度

朱宝才，王文龙，梁文俊

(1.山西农业大学 林学院，山西 太谷 030801； 2.西北农林科技大学 水土保持研究所黄土高原土壤侵蚀与旱地农业国家重点实验室，陕西 杨凌 712100； 3.中国科学院 水利部水土保持研究所，陕西 杨凌 712100)



用ArcGIS的Python脚本计算采煤沉陷区土地损毁程度

朱宝才1,2，王文龙2,3*，梁文俊1

(1.山西农业大学 林学院，山西 太谷 030801； 2.西北农林科技大学 水土保持研究所黄土高原土壤侵蚀与旱地农业国家重点实验室，陕西 杨凌 712100； 3.中国科学院 水利部水土保持研究所，陕西 杨凌 712100)

摘要：为计算井工煤矿采煤沉陷区土地损毁程度，利用ArcGIS的Python脚本语言编制软件，实现通过输入相同坐标系的下沉等值线、附加倾斜等值线及水平变形等值线3个DWG格式文件和土地利用现状SHP格式文件，直接输出土地损毁程度的SHP格式文件，从而获得沉陷区土地损毁程度预测图。通过对三元吉祥煤矿的实例分析，该软件亦可直接用于其他井工煤矿开采沉陷区土地损毁程度计算，且方便基层技术人员使用和推广。

关键词：采煤沉陷区；土地损毁；预测；GIS；Python

2011年5月，国土资源部制定《土地复垦方案编制规程》，共分7个部分，其中第3部分：井工煤矿(TD/T 1031.3 —2011)附录A 给出井工煤矿开采土地损毁预测方法，已有成熟软件[1～8]计算下沉、倾斜变形和水平变形；附录 B 给出采煤沉陷土地损毁程度分级参考标准，尚无成熟软件进行土地损毁程度分级计算，要在GIS软件中结合土地利用现状图手动分析，效率较低，且容易出错。为快速准确地计算采煤沉陷土地损毁程度，本研究以三元吉祥煤矿为例根据土地复垦方案编制规程的采煤沉陷土地损毁程度分级参考标准，基于井工煤矿开采土地损毁预测和土地利用现状数据，用ArcGIS 脚本语言Python编制“采煤沉陷土地损毁程度分级计算软件V1.0”(已登记计算机软件著作权，证书号：软著登字第1086259号)实现采煤沉陷区土地损毁程度计算。

1研究区概况与研究方法

1.1研究区概况

三元吉祥煤矿地理坐标为东经 113°03′09″～113°05′11″，北纬 36°22′07″～36°22′56″，海拔最高1 029.92 m，最低899.40 m，位于山西省长治市郊区黄碾镇，太行山中段西侧，长治盆地东部，总体北高南低，广为第四系黄土覆盖，北部有少许基岩出露，沟谷较发育，无常年性河流，雨季洪水汇入浊漳河。属于暖温带半湿润大陆性季风气候，四季分明，雨热同季，温和适中。据长治市气象局1961—2010年资料统计，年均气温 9.1 ℃，最高气温37.6 ℃，最低气温为-23.0 ℃。年均降雨量595 mm，最大降雨量832.9 mm(1971年)，最小降雨量340.9 mm(1965年)，雨季多集中在七、八2月，年均蒸发量1 558 mm，年平均无霜期170 d，冰冻期为11月至次年3月，最大冻土深度75 cm。一般冬季多西北风，夏季多东南风，最大风力10级。土壤类型为褐土性壤土，农作物主要是玉米和小麦，乔木主要是油松和刺槐，灌木以荆条、黄刺玫和胡枝子等为主，草本以白羊草草丛和蒿类草丛为主。

1.2数据来源

土地利用现状图数据源自第二次全国土地调查数据和2011年土地利用变更数据，为3度带高斯投影超图格式文件，在超图中转换文件格式为SHP格式，并重命名为LandUse_3.shp。

井工煤矿开采土地损毁预测数据采用概率积分法，通过输入与土地利用现状数据坐标系统一致的地形图、煤层底板等高线和工作面平面位置等数据，应用中国矿业大学“任意形多工作面多线段开采沉陷预计系统(MSPS)[1]”获得，包含10 mm下沉等值线和起点与间距均为0.5 m的下沉等值线文件命名为XiaChen.DWG。包含起点为0 mm·m-1，间距为2 mm·m-1的附加倾斜等值线文件命名为QingXie.DWG，包含起点为0 mm·m-1，间距为1 mm·m-1的水平变形等值线文件命名为ShuiPing.DWG。

以上数据皆存储于c:data中。

1.3数据处理

依据《土地复垦方案编制规程—第3部分：井工煤矿》推荐的“采煤沉陷土地损毁程度分级参考标准”，应用极限条件法确定各地类各等级土地损毁范围，再按原土地利用现状图、轻度以上损毁范围，中度以上损毁范围和重度损毁范围自下而上的顺序叠加即得土地损毁程度图，其中：

轻度以上损毁范围：依据国家煤炭工业局2000年制定的《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》确定10 mm下沉等值线包含区域。

水田中度以上损毁范围：{下沉>1.0 m区域}∪{水平变形>3.0 mm·m-1区域}∪{附加倾斜>4.0 mm·m-1区域}；水浇地中度以上损毁范围：{下沉>1.5 m区域}∪{水平变形>4.0 mm·m-1区域}∪{附加倾斜>6.0mm·m-1区域}；旱地中度以上损毁范围：{下沉>2.0 m区域}∪{水平变形>8.0 mm·m-1区域}∪{ 附加倾斜20.0 mm·m-1区域}；林地、草地中度以上损毁范围：{下沉>2.0 m区域}∪{水平变形>8.0 mm·m-1区域}∪{ 附加倾斜>20.0 mm·m-1区域}。

水田重度损毁范围：{下沉>2.0 m区域}∪{水平变形>6.0 mm·m-1区域}∪{附加倾斜>10.0 mm·m-1区域}；水浇地中度损毁范围：{下沉>3.0 m区域}∪{水平变形>8.0 mm·m-1区域}∪{ 附加倾斜>12.0 mm·m-1区域}；旱地中度损毁范围：{下沉>5.0 m区域}∪{水平变形>16.0 mm·m-1区域}∪{附加倾斜>40.0 mm·m-1区域}；林地、草地中度损毁范围：{下沉>6.0 m区域}∪{水平变形>20.0 mm·m-1区域}∪{附加倾斜>50.0 mm·m-1区域}

1.4Python实现

Python是开源编程语言，ArcGIS9.0开始引入Python。封装大部分ArcGIS 空间分析工具于Arcgisscripting 站点包中，从ArcGIS10.0 开始，站点包更名为Arcpy，Python脚本可独立于ArcGIS运行[9～14]。

1.导入ArcPy。

ArcPy提供所有地理处理工具使用 Python 语言操作的接口，并为处理和查询 GIS 数据提供函数和类，如env 类包含所有地理处理环境。通过导入env类，设置工作空间为c:/data文件夹，代码为：

import arcpy

from arcpy import env

env.workspace = “C:/data”

2.创建个人地理数据库存放临时文件

在C:/data创建名为Temp.gdb的个人地理数据库存放地理处理的临时文件，方便用后统一删除，代码为：

arcpy.CreateFileGDB_management(env.workspace, “Temp.gdb”)

3.分别转换CAD格式的下沉等值线文件XiaChen.DWG，水平变形等值线文件ShuiPing.DWG和附加倾斜等值线文件QingXie.DWG为地理数据库Temp.gdb的要素类XiaChen，ShuiPing和QingXie。以转换CAD格式下沉等值线为地理数据库为例，代码如下：

arcpy.CADToGeodatabase_conversion(“XiaChen.DWG”, “Temp.gdb”, “XiaChen”, “1000”)

拷贝土地利用现状LandUse_3.shp为Temp.gdb的要素类LandUse_3。代码为：

arcpy.CopyFeatures_management(“LandUse_3.shp”,“Temp.gdb/LandUse_3”)

4.用要素类XiaChen，ShuiPing和QingXie等值线创建临时图层，拟存放临界下沉等值线(值为10, 1000, 1500, 2000, 3000, 5000和6000)、临界水平变形等值线(值为3, 4, 6, 8, 16, 20)和临界附加倾斜等值线(值为4,6,10,12,20,40,50)，图层名命名采取“要素类名”+“_”(下划线)+“临界值”的方法，拟存放下沉值为10的等值线图层命名为“XiaChen_10”；拟存放水平变形值为3的图层命名为“ShuiPing_3”；拟存放附加倾斜值为4的图层命名为“QingXie_4”。其代码示例如下：

arcpy.MakeFeatureLayer_management(“Temp.gdb/XiaChen/Polyline”, “XiaChen_10”)

arcpy.MakeFeatureLayer_management(“Temp.gdb/ShuiPing/Polyline1”, “ShuiPing_3”)

arcpy.MakeFeatureLayer_management(“Temp.gdb/QingXie/Polyline2”, “QingXie_4”)

用要素类LandUse_3创建临时图层，拟存放水田、水浇地、旱地，林地和草地及全部地类的图层分别命名为LandUse_011、LandUse_012、LandUse_013和LandUse_03_04和LandUse。

5.通过属性“Elevation”选择临界下沉等值线、临界水平变形等值线和临界附加倾斜等值线，图层名中包含“临界值”即为临界等值线值，XiaChen_10选择的属性即为“Elevation”=10。因水平变形和附加倾斜值用正负号表示变形方向，故对其等值线值取绝对值。其代码示例如下：

arcpy.SelectLayerByAttribute_management(“XiaChen_10”, “NEW_SELECTION”, ‘“Elevation”= 10’)

arcpy.SelectLayerByAttribute_management(“ShuiPing_3”, “NEW_SELECTION”, ‘Abs(“Elevation”)=3’)

arcpy.SelectLayerByAttribute_management(“QingXie_4”, “NEW_SELECTION”, ‘Abs(“Elevation”) = 4’)

通过地类编码属性“DLBM”选择地类，011为水田，012为水浇地，013为旱地，031为有林地，以上各类通过新选择，再将032、033、041、042和043添加选择到LandUse_03_04图层中。其代码示例如下：

arcpy.SelectLayerByAttribute_management(“LandUse_03_04”, “NEW_SELECTION”, ‘CAST(“DLBM” AS INTEGER) = 031’)

arcpy.SelectLayerByAttribute_management(“LandUse_03_04”, “ADD_TO_SELECTION”, ‘CAST(“DLBM” AS INTEGER) = 032’)

6.将临界等值线转换为临界面，输入文件为已通过对应属性选择的临界等值线图层，输出文件命名为“polygon_”+原临界等值线文件名，如XiaChen_10等值线转换为面后，命名为“polygon_XiaChen_10”，其代码示例如下：

arcpy.FeatureToPolygon_management(“XiaChen_10”, “Temp.gdb/polygon_XiaChen_10”, “”, “NO_ATTRIBUTES”, “”)

7.确定轻度及轻度以上损毁范围。在polygon_XiaChen_10属性表中添加损毁等级字段“SHDJ”，字段类型为SHORT，在字段中可存储9位。并将“SHDJ”字段统一赋值为“1”，即轻度损毁。代码实现如下：

arcpy.AddField_management(“Temp.gdb/polygon_XiaChen_10”, “SHDJ”, “SHORT”, 9, “”, “”,“SHDJ”, “NULLABLE”)

arcpy.CalculateField_management(“Temp.gdb/polygon_XiaChen_10”, “SHDJ”, “1”, “”)

8.分地类确定中度及中度以上和重度可能损毁范围。依据下沉、水平变形、附加倾斜临界值联合确定水田、水浇地、旱地、林草地中度及中度以上可能损毁范围，要素类分别命名为Land011_Class_2，Land012_Class_2，Land013_Class_2和Land03_04_Class_2。水田、水浇地、旱地、林草地重度可能损毁范围要素类分别命名为Land011_Class_3，Land012_Class_3， Land013_Class_3和Land03_04_Class_3。其实现代码示例如下：

arcpy.Union_analysis ([“Temp.gdb/polygon_XiaChen_1000”,“Temp.gdb/polygon_ShuiPing_3”, “Temp.gdb/polygon_QingXie_4”], “Temp.gdb/Land011_Class_2”,“ONLY_FID”,“”)

在以上各要素属性表中添加损毁等级字段“SHDJ”，设置同前，并将损毁等级为中度及中度以上要素(即名称中有Class_2)的“SHDJ”字段统一赋值为“2”，损毁等级为重度要素(即名称中有Class_3)的“SHDJ”字段统一赋值为“3”。

9.融合，为使制图中同一地块不被分割成具有相同属性的多个地块，将以上要素按“SHDJ”字段进行融合，融合后的要素名在原要素名上加“_dissolve”。其代码示例如下：

arcpy.Dissolve_management(“Temp.gdb/Land011_Class_2”, “Temp.gdb/Land011_Class_2_dissolve”, “SHDJ”, “”, “”, “”)

10.分地类确定中度和重度可能损毁范围。更新要素类属性，用水田重度可能损毁范围Land011_Class_3_dissolve更新中度及中度以上可能损毁范围Land011_Class_2_dissolve，得到水田中度和重度可能损毁范围Land011_Class_2_3，代码实现如下：

arcpy.Update_analysis(“Temp.gdb/ Land011_Class_2_dissolve”, “Temp.gdb/Land011_Class_3_dissolve”,“Temp.gdb/Land011_Class_2_3”)

同样方法，分别得到水浇地、旱地、林草地的中度和重度可能损毁范围要素类Land012_Class_2_3，Land013_Class_2_3和Land03_04_Class_2_3。

11.标识(Identity)要素类，用轻度及轻度以上损毁范围要素Temp.gdb/polygon_XiaChen_10标识土地利用Temp.gdb/LandUse，得到轻度及轻度以上损毁土地Temp.gdb/LandUse_Damage_1，实现代码如下：

arcpy.Identity_analysis (“Temp.gdb/LandUse”, “Temp.gdb/polygon_XiaChen_10”, “Temp.gdb/LandUse_Damage_1”)

用水田中度和重度可能损毁范围Land011_Class_2_3标识水田土地利用LandUse_011，得到水田中度和重度损毁范围LandUse011_Damage_2_3。实现代码如下：

arcpy.Identity_analysis (“LandUse_011”, “Temp.gdb/Land011_Class_2_3”, “Temp.gdb/LandUse011_Damage_2_3”)

同样方法分别得到水浇地、旱地、林草地的中度和重度土地损毁要素类Land012_ Damage _2_3，Land013_ Damage _2_3和Land03_04_ Damage _2_3。

12.逐地类更新(Update)，用水田的中度和重度土地损毁要素类LandUse011_Damage_2_3更新轻度及轻度以上损毁土地LandUse_Damage_1，再用水浇地的中度和重度土地损毁要素类Land012_ Damage _2_3更新其结果，依次Land013_ Damage _2_3和Land03_04_ Damage _2_3更新前一更新结果，实现代码示例如下：

arcpy.Update_analysis (“Temp.gdb/LandUse_Damage_1”, “Temp.gdb/LandUse011_Damage_2_3”, “Temp.gdb/LandUse_Damage_011OK”)

arcpy.Update_analysis (“Temp.gdb/LandUse_Damage_013OK ”, “Temp.gdb/LandUse03_04_Damage_2_3”, “LandUse_Damage”)

将上述代码保存为脚本文件Damage.py。

1.5软件运行

在ArcGIS10.2可运行的状态下，用Python2.7 打开本脚本程序文件Damage.py，按F5键运行脚本。

为获得数据处理所用时间，在开始运行软件时获取系统时间，单位是秒，数据处理结束后再次获取系统时间，两者之差即为数据处理总时间。本次数据处理在Intel(R) Core(TM) i5-4210Y CPU@1.50 GHz四核处理器，4.00 GB内存，128 GB固态硬盘的条件下运行240 s。

计算结果土地损毁程度文件LandUse_Damage.shp存储于C:data文件夹中，中间过程数据存储于Temp.gdb地理数据库中，LandUse_Damage.shp文件相对于LandUse_3.shp增加了三个字段，分别是SHDJ(损毁等级)、Shape_Leng(周长)和Shape_Area(面积)，SHDJ字段所列数字即为该小班采煤沉陷土地损毁程度，其中0代表无损毁；1代表轻度损毁；2代表中度损毁；3代表重度损毁。Shape_Area字段为该小班面积。

2结果与分析

基于“采煤沉陷土地损毁程度分级计算软件V1.0”计算的三元吉祥煤矿土地损毁程度预测面积统计见表1，由于篇幅所限，未附土地损毁程度预测图。各地类轻度损毁总面积为361.47 hm2，占损毁总面积的63.20%，重度损毁总面积123.42 hm2，占损毁总面积的21.58%，可见三元吉祥煤矿开采沉陷区土地损毁以轻度为主。损毁的耕地全部为旱地，面积为257.88 hm2，占总损毁土地面积的45.09%，其中轻度损毁占23.85%，中度损毁占7.80%，重度损毁占13.44%。损毁林地包括有林地和其他林地，面积为121.74 hm2，占总损毁土地面积的21.28%，其中轻度损毁占11.79，中度损毁占5.17%，重度损毁占4.33%。损毁草地全部为其他草地，面积为66.32 hm2，占总损毁土地面积的17.67%，其中轻度损毁占11.60%，中度损毁占2.26%，重度损毁占3.81%。耕地林地草地共损毁480.66 hm2，占损毁总面积的84.04%，主要原因是矿区内以耕地林地草地为主。

表1三元吉祥煤矿采煤沉陷区土地损毁程度预测面积统计表/hm2

Table 1Sanyuan Jixiang coal mining subsidence land damage statistical prediction of surface area

地类名称Landuse轻度损毁Lightdamage中度损毁Mediumdamage重度损毁Severedamage总计Total旱地136.4344.6076.85257.88有林地31.2117.4912.5161.21其他林地36.2012.0812.2560.53其他草地66.3212.9021.82101.04铁路用地3.140.000.003.14农村道路0.300.000.000.30设施农用地1.260.000.001.26沼泽地0.710.000.000.71裸地6.050.000.006.05城市1.420.000.001.42村庄47.260.000.0047.26采矿用地30.830.000.0030.83风景名胜及特殊用地0.340.000.000.34总计361.4787.07123.42571.96

基于“采煤沉陷土地损毁程度分级计算软件V1.0”计算的三元吉祥煤矿土地损毁程度预测小班731个，基于手动ArcGIS预测方法，获得土地损毁程度小班也是731个，将两种方法获得的土地损毁程度小班分别按面积从大到小排序，逐一对比，发现对应地类完全相同，对应面积差小于0.001 m2，相对误差小于0.1%。

3结论和讨论

基于“采煤沉陷土地损毁程度分级计算软件V1.0”进行煤矿土地损毁程度预测可达到手动ArcGIS分析相同的精度。利用“采煤沉陷土地损毁程度分级计算软件V1.0”能快速地计算采煤沉陷区土地损毁程度，可直接用于其他井工煤矿开采沉陷区土地损毁程度计算，为土地复垦方案，生态环境恢复治理方案及水土保持方案等提供准确的土地损毁程度数据。

“采煤沉陷土地损毁程度分级计算软件V1.0” 处理时间大量缩短，技术要求显著降低，可不打开ArcGIS，直接在Python2.7中运行，运行前只需把土地利用现状、下沉等值线、附加倾斜等值线和水平变形等值线文件按指定格式和名称存储于指定位置即可，方便简单，降低了采煤沉陷区土地损毁程度计算门槛，方便基层技术人员使用和推广。

参考文献

[1]蔡怀恩. 开采沉陷预计的方法及发展趋势[J]. 露天采矿技术, 2007, 22(4):43-44.

[2]康建荣, 王金庄, 温泽民. 任意形多工作面多线段开采沉陷预计系统(MSPS)[J]. 矿山测量, 2000, 28(1):24-27.

[3]李明, 李琰庆, 王红梅. 任意形状工作面开采沉陷预计系统开发[J]. 矿业安全与环保, 2008, 35(5):18-21.

[4]李培现, 谭志祥, 齐公玉, 等. 基于MATLAB的开采沉陷预计系统[J]. 中国矿业, 2008,17(11):72-76.

[5]刘虎, 吴侃, 陈冉丽. 基于AutoCAD的开采沉陷预计分析软件研究[J]. 矿业安全与环保, 2009, 36(2):22-23.

[6]谢和平, 周宏伟, 王金安, 等. FLAC在煤矿开采沉陷预测中的应用及对比分析[J]. 岩石力学与工程学报, 1999, 18(4):29-33.

[7]张庆松, 高延法, 刘松玉, 等. 基于粗集与神经网络相结合的岩移影响因素分析与开采沉陷预计方法研究[J]. 煤炭学报, 2004, 29(1):22-25.

[8]柴华彬, 邹友峰, 袁占良, 等. SuperMap系统在开采沉陷预计分析中的应用[J]. 西安科技大学学报, 2005,25(1):52-56.

[9]方圣辉, 张玉贤, 佃袁勇, 等. 基于Python的ArcGIS地理数据批处理[J]. 测绘与空间地理信息, 2015, 38(1):1-2.

[10]焦洋, 邓鑫, 李胜才. 基于Python的ArcGIS空间数据格式批处理转换工具开发[J]. 现代测绘, 2012, 35(3):54-55.

[11]李强, 白建荣, 李振林, 等. 基于Python的数据批处理技术探讨及实现[J]. 地理空间信息, 2015, 13(2):54-56.

[12]彭海波, 向洪普. 基于Python的空间数据批量处理方法[J]. 测绘与空间地理信息, 2011, 34(4):81-82.

[13]邵保华, 田学志. 谈Python在Arcgis地理处理中的应用[J]. 林业勘查设计, 2012, 41(2):99-100.

[14]朱海金, 翁伟琳. 基于python数字高程模型地形数据批量提取[J]. 地理空间信息, 2013, 11(5):136-138.

(编辑：马荣博)

Calculate land damage degree in the coal mining subsidence area by Python, ArcGIS script language

Zhu Baocai1, 2, Wang Wenlong2, 3*， Liang Wenjun1

(1.CollegeForestry,ShanxiAgriculturalUniversity,Taigu030801,China; 2.StateKeyLaboratoryofSoilandErosionandDrylandFarmingontheLoessPlateau,InstituteofSoilandWaterConservation,NorthwestA&FUniversity,Yangling712100,China; 3.InstituteofSoilandWaterConservation,ChineseAcademyofSciencesandMinistryofWaterResources,Yangling712100,China)

Abstract:For calculating land damage degree in the coal mining subsidence area, the Python that attached to ArcGIS, a script language compilation software, were used to predict the land damage degree area of every land use type. Land damaged degree result characterized as the SHP file format can be outputted directly and then the predicted figures of land damaged can be obtained, through the way that the subsidence contour, additional tilting contour and horizontal deformation contour files with DWG format and the land use file with SHP format were inputted altogether in the Python software. The satisfactory result was obtained from SanYuanJiXiang Coal Mine. Therefore, the software can be applied directly to calculate land damaged degree in other coal mining subsidence area. It is promoted easily and convenient for the application of basic technician.

Key words:Coal mining subsidence area; Land damage; Prediction; GIS; Python

中图分类号：F301.24

文献标识码：A

文章编号：1671-8151(2016)05-0371-06

基金项目：国家自然科学基金项目(31500523，41501201，41571275)，中国科学院西部行动计划项目(KZCX2-XB3-13-4)

作者简介：朱宝才(1978-)，男(汉)，内蒙古巴林左旗人，讲师，博士研究生，研究方向：土壤侵蚀与水土保持*通讯作者：王文龙，博士，研究员，博士生导师。Tel：13669222856；E-mail：wlwang@nwsuaf.edu.cn

收稿日期：2015-12-18 修回日期：2016-01-28