采煤塌陷积水区面积无人机采集方法

迟臣鑫 陈伟强 朱鹏程 查剑锋 殷和健 马晨阳

（1.中国矿业大学环境与测绘学院，江苏 徐州 221116；2.江苏省资源环境信息工程重点实验室，江苏 徐州 221116；3.皖北煤电集团有限责任公司五沟煤矿，安徽 淮北 235131）

雨水囤积与地下水汇流会在采煤塌陷区形成积水，会破坏矿区及附近生态环境［1-2］、造成耕地、道路或建筑物的破坏。据统计，截至2018年，我国两淮地区的采煤塌陷区面积约300 km²，积水面积占塌陷区面积的30%～50%［3］。近年来，采煤塌陷积水区的生态修复得到了学术界的普遍关注［4-6］，精确获取面积这一基础数据，对于有效开展采煤塌陷积水区生态修复［7］等工作具有重要意义。

采煤塌陷积水区面积获取方法通常有RTK实测法和遥感法，后者最为常用。范忻等［8］提出了一种改进的P-WSVM方法，用于计算塌陷积水区面积，该方法受到遥感数据分辨率和波段的影响，提取精度为87.76%；王安妮［9］利用NCIWI指数法，以Landsat影像为数据源提取了2000—2018年的沛北矿区积水区面积；胡振琪等［10］采用基于小波变换的多尺度边缘检测算法提取了塌陷积水区边界。一般而言，采煤塌陷积水区具有面积现势性强、地形复杂及测区形状不规则的特点，传统方法存在的不足有：①RTK实测法适用于小区域的采煤塌陷积水区面积计算，若要保证面积计算精度，则需提高RTK实地采点间隔，工作量会随之增加；②考虑丰水期、枯水期时采煤塌陷积水区面积有差异，并受到遥感影像分辨率的影响［11］，遥感法适用于大范围的采煤塌陷积水区的面积计算，且高分辨率遥感影像获取成本较高，计算范围较小的采煤塌陷积水区面积会产生数据冗余和经费浪费。为此，本研究提出了一种以无人机影像为数据源，利用耦合滤波方法消除噪声、自动提取水域轮廓并最终获取采煤塌陷积水区面积的方法。

1 塌陷积水区面积信息无人机采集方法

无人机采集地表信息的一般流程有航线规划、布设像控点和提取信息。由于塌陷积水区随煤炭开采逐步变化，其面积在不同季节会有差异并且水域边界不规则，故本研究提出了如图1所示的面积采集思路。

具体实施步骤为：①围绕采煤塌陷积水区布设像控点并用无人机拍摄正射影像［12］，保证航向、旁向重叠率在70%以上；②对获得的正射影像进行预处理，包括灰度化、去噪和二值化处理等；③提取并处理正射影像上积水区的轮廓，将轮廓点位坐标转换为实地坐标，采用向量叉乘法计算积水区面积；④对采煤塌陷积水区面积的计算精度进行评估。

2 塌陷积水区面积信息提取的关键技术

2.1 影像数据处理

2.1.1 灰度化

灰度化是影像二值化处理的前提，为保留采煤塌陷积水区影像灰度化后的对比度，采用文献［13］提出的基于最大加权投影求解的灰度化算法进行影像灰度化处理。

2.1.2 噪声处理

常见的图像噪声有高斯噪声和椒盐噪声等。高斯噪声是指概率密度函数服从高斯分布的噪声，椒盐噪声是指图像传感器、传输信道及解码处理在影像获取、传输等过程中产生的黑白相间的亮暗点噪声［14］。图像去噪的经典方法有中值滤波、高斯滤波、均值滤波、高通滤波等。高斯滤波适用于处理高斯噪声［15］；均值滤波同样对高斯噪声的去除效果较好，但会损失图像的边缘信息［16］；高通滤波可以抑制低频噪声，对图像进行锐化；中值滤波较适合去除椒盐噪声。目前，已有较多学者对上述滤波方法进行了改进［17-20］，无人机采集的采煤塌陷积水区影像中存在多种噪声，单一滤波方法降噪效果有限，采用滤波耦合去噪具有较好的效果。本研究耦合去噪采用中值滤波、高通滤波算法实现。其中，中值滤波选择5×5的像素模板，是因为采用较小尺寸模板（3×3）虽然能大量保留图像细节信息，但图像平滑效果较差；采用较大尺寸模板（7×7）虽然平滑效果好，但会丢失图像细节信息；高通滤波算法通过sobel算子实现，可达到突出轮廓的目的。图2为采用不同滤波方法得到的采煤塌陷积水区轮廓。

分析图2可知：采用本研究耦合滤波去噪思路处理后，图像中噪声得到了有效抑制，水体轮廓信息提取较为精准，效果优于单纯采用高通滤波算法、中值滤波算法以及不进行任何滤波处理情况下的效果。图2中，除了图2（b）中噪声点较多，无法计算轮廓面积外，其余3类算法处理后的轮廓面积分别为674.34 km2、779.10 km2、726.39 km2，只有滤波耦合去噪算法得到的轮廓面积与依据矿区实测数据的Pix4Dmapper软件计算值（716.82 km2）最为接近，可见，本研究提出的滤波算法有助于减少图像噪声，突出积水区的轮廓信息。

2.1.3 二值化处理

二值化处理可以达到减少数据量并突出采煤塌陷积水区轮廓信息的目的，常用的阈值选取方法有［21］：三角法求阈值［22］、大津法［23］、平均值阈值法以及固定阈值法。采用该类方法生成的二值化图如图3所示。

由图3可知：采用三角法确定阈值时，水体轮廓清晰，边缘平滑且无明显噪声点，与实际水域的轮廓最为接近；采用大津法确定阈值时，轮廓边缘噪声点较多，无法准确区分水域边缘的杂草、湿地等干扰信息，提取的水域轮廓面积较实际偏大；采用平均值法选取阈值时，轮廓边缘不平滑，有较多的噪声点被提取出；采用固定阈值法提取的轮廓边缘较为平滑，轮廓提取速度较快，噪声点也较少，提取的轮廓面积与实际面积相接近，但对于不同的图像，较难选取一个通用的固定阈值。故本研究采用三角法进行阈值选取。

2.2 塌陷积水区面积计算方法

2.2.1 提取水域轮廓及坐标

OpenCV中的Findcontours函数具有提取影像轮廓的功能，将处理后的二值化正射影像输入，则会生成完整的嵌套轮廓层次结构并返回Contours列表。

2.2.2 坐标转换

轮廓坐标为图上坐标，要计算轮廓面积，需将图上坐标转换为实地坐标。图上坐标转换为实地坐标的公式为

式中，trans(0)、trans(3)代表图像左上角点的实地坐标；trans(1)代表像元宽度；trans(5)代表像元高度；若图像指北，则参数trans(2)、trans(4)默认为0；m、n为水域轮廓上各个点的像素坐标；x、y为水域轮廓上各个点的实地坐标。

2.2.3 面积计算

三角形的面积Area(i)等于其两边向量积（叉乘）的模的1/2，本研究将任意n边形分解为n-2个三角形并求和，来计算不规则多边形的面积Y，计算公式为

式中，Area(i)为第i个三角形的面积；n为轮廓的顶点数（边数）。

读取水域轮廓上各个点的x、y坐标，并将其转换为顶点Zp(p=1,2,3,…,n)，即水体的轮廓是顶点为Zp的任意多边形，顶点坐标按逆时针方向排列，分别为(x1,y1),(x2,y2),…,(xn,yn)。由于提取的坐标点数量较多，故采煤塌陷积水区可近似为不规则多边形，面积计算公式为

式中，x,y为轮廓点的坐标。

2.3 塌陷积水区面积计算精度评估方法

塌陷积水区面积计算精度受到无人机正射影像获取精度、正射影像处理效果、水域轮廓坐标点提取精度和面积计算精度等因素影响。因水域轮廓坐标点较多，故其精度对面积计算的影响最大。本研究依据轮廓坐标点的精度来评估面积计算精度。实现步骤有：①生成N(N>100)组轮廓坐标点的高斯随机数，赋予轮廓坐标点不同的方差；②计算N组带有随机方差的轮廓面积；③比较每组面积的期望值与依据相关算法计算值的差异。不同坐标点位测量误差对应的面积误差如图4所示。

3 实例分析

3.1 塌陷积水区概况

采煤塌陷积水区位于某矿开采区，如图5所示，采煤引起地表塌陷，雨水囤积与地下水汇流导致该区域形成了3处采煤塌陷积水区。3处塌陷积水区内部均无水草覆盖，水体边缘有少量水草和树木，水域1、2轮廓形状不规则，用作农田灌溉蓄水池，水域3为连通复杂轮廓，用作鱼塘养殖。为满足该矿对塌陷积水区生态修复及安全防护的要求，应用所提方法对3处积水区的面积进行了计算分析。

3.2 数据采集

在以3处积水区为中心的矩形区域上分别布设了8个像控点，采用大疆Phantom4 RTK无人机获取塌陷积水区的正射影像，用Pix4Dmapper及Global Mapper软件提取并处理正射影像信息。

3.3 塌陷积水区面积提取

根据本研究方法实现思路，使用python语言编写程序，依据无人机影像自动计算采煤塌陷积水区的面积。以水域2为例，该处水域轮廓如图6所示，水域2轮廓点位提取结果如图7所示。

3.4 塌陷积水区面积测量精度分析

本研究分别采用Pix4Dmapper软件、RTK实测法以及本研究方法计算了塌陷积水区的面积，结果如表1所示。

由表1可知：本研究所提方法与其他两种方法计算的水域面积相对误差在6%以内；所提方法计算的水域1面积大于Pix4Dmapper软件，该方法计算的水域2、水域3的面积均小于RTK实测法和Pix4Dmap⁃per软件，原因在于水域2、水域3的水体边缘有较多树木、杂草等干扰，水域1的水体边缘树木、杂草等干扰较少。所提方法在提取水体轮廓时会将水体边缘的树木、杂草等干扰信息排除，故计算的面积偏小；当水体边缘无树木、杂草等干扰信息时，轮廓提取较为精细，故计算的面积偏大。

4 结论

（1）为满足采煤塌陷积水区生态修复需求，提出了一种基于无人机影像的采煤塌陷积水区面积采集方法。结合某矿工程实例，对该方法的适用性进行了分析。结果表明：该方法与商用化软件Pix4Dmapper及RTK实测法相比，采煤塌陷积水区面积计算误差在6%以内，表明该方法对于采煤塌陷积水区面积计算有一定的适用性，可为矿山生态修复等工作提供可靠依据。

（2）由于无人机航摄存在续航时间较短、像控点布设数量与测区面积成正比的不足，故该方法适用于中小面积的采煤塌陷积水区面积计算。相对于RTK实测法和遥感法，该方法工作量较少，且时效性较强。