基于无人机的高潜水位煤矿区沉陷耕地提取方法比较

胡 晓，李新举

(1.山东农业大学 资源与环境学院，山东 泰安 271018; 2.山东农业大学 信息科学与工程学院，山东 泰安 271018)

高潜水位煤矿区主要分布于我国中东部平原的河北省、河南省、山东省、安徽省和江苏省5省[1-2]，且煤粮复合面积大[3-4]。采煤沉陷后地下水位上升到地表标高以上，导致大量高产优质耕地常年积水或季节性积水，造成耕地大幅度减产甚至绝产[5-6]，严重影响了煤矿区农民的生产生活。因此准确、快速提取沉陷耕地的信息，不仅可以为耕地损害范围的测算提供科学依据，而且关系到矿区居民的切实利益，是维护矿区社会稳定的需要。

目前，信息的提取大多依靠卫星遥感技术手段，一些学者利用卫星遥感技术开展了煤矿区积水沉陷耕地的提取研究[7-9]，但是卫星遥感影像分辨率低、获取周期长、易受云雾遮挡[10]，在非积水沉陷耕地信息提取方面则显得捉襟见肘。近几年飞速发展的无人机遥感技术为解决这个问题提供了新的思路，作为卫星遥感的有益补充，无人机具有响应快、周期短、易操作、成本低等特点[11]，传感器获得的高分辨率的遥感影像在表达地物几何、纹理、拓扑关系等特征参量方面更加细致[12]，增强了对地物类型的识别能力，已经在耕地作物信息提取[13-14]等方面得到较为广泛的应用。已有研究中，无人机搭载的多为机载高光谱和多光谱传感器，波段包含可见光-近红外波段，昂贵的传感器制约了无人机遥感的发展和应用[15]。

仅含可见光波段的无人机影像获取简便、成本低廉，实际使用也更为广泛[16]。基于无人机可见光影像的信息提取方法大致有两类。一是采用图像增强的方法，采用各种增强算法，可以提高图像的目视解译效果，方便人工目视解译、图像分类中的样本选取等，以便使结果图像比原始图像更适合于特定的应用要求。如李宗南等[17]利用小型无人机获取的可见光图像研究了灌浆期玉米倒伏的图像特征和面积提取方法，分别利用基于色彩特征和评选出的纹理特征两种图像增强方法提取倒伏玉米面积，结果显示，基于红、绿、蓝色均值纹理特征提取倒伏玉米面积的误差显著低于基于色彩特征提取方法。二是采用基于可见光三波段的植被指数提取的方法，如汪小钦等[18]借鉴归一化植被指数NDVI 的构造原理及形式，提出了一种综合利用红、绿、蓝3个可见光波段的归一化植被指数-可见光波段差异植被指数VDVI(Visible-band difference vegetation index)，并与其它一些可见光植被指数进行了比较，结果表明，VDVI 指数提取精度高于其他可见光波段植被指数，且阈值在0 附近，较易确定。VDVI是一种基于波段运算的植被指数，除此之外，可见光植被指数还包括基于单一波段和基于可见光植被指数运算两种形式。

以上相关研究大都以健康农作物为对象，以高潜水位煤矿区沉陷耕地农作物作为对象的研究较少，方法主要以图像增强或者植被指数单一提取为主，并未对这两种方法进行比较与分析。利用无人机的可见光波段，以高潜水位煤矿区作为研究区，开展沉陷耕地信息的提取研究还处于初级阶段，而且方法的有效性和实用性还有待进一步验证。因此本研究以仅含可见光波段的无人机影像为数据源，在前人的研究基础上，对高潜水位煤矿区沉陷耕地提取方法开展研究比较，以期能够获得一种准确、有效、低成本的沉陷耕地信息提取方法，扩大无人机遥感在高潜水位煤矿区应用的深度和广度，为沉陷耕地的快速准确识别、动态变化监测及土地复垦等提供数据基础和技术支撑。

1 材料与方法

1.1 试验区概况及数据获取

试验区隶属于兴隆庄煤矿，位于山东省济宁市兖州区兴隆庄镇境内，其井田位于兖州煤田中北部，京沪铁路南北纵贯井田，是典型的高潜水位煤矿区。影像获取时间为2018-04-13T9:00—11:00,天气晴朗无风。采用的是科卫泰X6L六旋翼无人机，飞行高度150 m，搭载红鹏相机获取试验区航片，利用Global Mapper数据处理软件对获取的航片进行校正、拼接处理，得到了试验区正射影像(图1)，影像包含红、绿、蓝3个波段，分辨率为0.03 m。

图1 试验区正射影像

试验区高程在37.8～54.3 m，整体地势从北向南逐渐降低，从图1可以看出，试验区南部及西南部已形成了沉陷积水区域，积水区在影像上呈现绿色和黑色，经实地调查，绿色水面为常年积水区，沉陷程度大、面积广;土壤背景颜色导致轻度沉陷积水区在影像上呈现黑色。耕地种植作物为冬小麦，积水区之间地表沉陷严重，沉陷积水引发的土壤养分流失和盐渍化，已经导致积水区周围耕地冬小麦生长严重滞后，密度大幅度降低，影像上呈现稀疏的浅绿色，与田块状耕地颜色反差明显，但是纹理仍然可见。本文以沉陷耕地为研究对象，从试验区中截取一块耕地区域作为研究区(图 1(a))。

1.2 无人机影像特征分析

为了准确提取研究区影像的耕地信息，分析无人机影像上冬小麦与非植被的光谱特性，通过观察它们在各波段间反射率的差异，构造植被指数，提取耕地信息。为了使典型地物光谱信息的统计更具代表性，对于每种地物分别选取10 个代表区域，统计典型地物在红、绿、蓝3 个波段中像元的统计特征值(包括最大值、最小值、平均值及标准偏差等)。选取平均值分析各地物在3 个波段的像元值总体差异，并引入标准差作为每个波段像元值的波动范围(式(1))，统计结果见表1。

(1)

表1 典型地物在蓝、绿、红波段的像元值差异

Table 1 Difference pixels value in blue,green,red band of typical land cover

典型地物蓝光波段像元值均值标准差绿光波段像元值均值标准差红光波段像元值均值标准差耕地96.35836.992109.38723.58486.48124.186道路216.59630.151203.12930.847188.32134.155积水区94.76524.82991.57422.27690.08720.691裸地165.94615.727154.96719.644147.80123.167芦苇202.58113.153185.81013.760164.95814.596

从表1中可以看出:① 对于耕地有绿光波段像元值>蓝光波段像元值>红光波段像元值，即绿光反射率最大，其次是蓝光波段与红光波段，这与健康绿色植被的光谱曲线不相吻合，这恰恰说明了采煤扰动下的沉陷导致潜水位相对升高，产生渍害影响冬小麦的正常生长，而且影响的范围较大;② 其它地物信息满足蓝光波段像元值>绿光波段像元值>红光波段像元值;③ 道路、裸地、芦苇3种地物在3个波段均无重合，这说明采用任一波段进行区分都较为容易实现;④ 芦苇不属于绿色植被，因此反射特征与绿色植被不相同;⑤ 耕地与积水区在3个波段均有一定重合，这是由潜水位上升形成的积水耕地造成的，绿光波段对冬小麦的反射率大于积水区，可以作为耕地提取的依据。

综上分析可以看出，由于耕地在3个波段上与积水区均存在重叠，因此，一方面，可以采用滤波方式进一步凸显耕地规则的纹理信息，或者采用色彩空间变换进一步增大各地物之间的色彩反差，滤波分析、色彩空间转换等处理有增强图像的作用，对耕地的提取有帮助;另一方面，可以通过绿光通道与红光、蓝光通道间的运算，或者通过植被指数间的进一步运算，增强耕地与周围地物的差异，实现耕地的准确提取。

1.3 色彩和纹理特征的统计与分析

在影像中建立典型地物的感兴趣区，计算影像的色彩特征和纹理特征并进行统计和分析，筛选出能够区分不同地物类型的关键特征参数，并以此为依据进行耕地信息的提取。

1.3.1色彩和纹理特征的统计与分析

为了进一步增强耕地与其它地物尤其是积水区的区别，本研究利用ENVI 5.3软件对获取的影像首先进行了HLS(hue，lightness，saturation)色彩空间转换，生成了一个增强的彩色合成图像，统计了色度、亮度及饱和度3项色彩特征，其次计算灰度共生矩阵，利用二阶矩阵纹理滤波得到了3种颜色的均值、方差、协同性、对比度、相异性、信息熵、二阶矩、相关性共24项纹理特征。其中，二阶概率滤波的窗口设置为7×7，空间相关性矩阵X和Y的变换分量为1和1，灰度质量级别为64。

为了使耕地的色彩和纹理信息统计更具代表性，从整个研究区中较为均匀的选取了20个耕地样本，15个积水区样本，利用ENVI 5.3对各项结果进行统计，最终得到了27项特征指标的均值和标准差，根据均值和标准差计算变异系数，再根据样本特征均值计算耕地和积水区之间的相对差异，结果见表2，具体计算方法：

恒山景区在旅游基础设施建设方面基础设施不到位，首先表现在停车场建设滞后、自驾车营地和自助游驿站等建设迟缓、标识导览系统短缺、游客服务中心功能不全面、厕所革命相对滞后等；二、景区景点的通信基础设施建设不到位，并没有覆盖景区全范围的主要通用网络和“旅游大数据云平台建设行动计划”，没有实现旅游服务和管理全程覆盖、全域覆盖。三、城市与景区之间、景区与景区之间的连接公路以及交通工具单一，交通通达能力和游客接待能力较弱。

(2)

(3)

式中，S为标准差;M为样本均值;M1为耕地区的像元均值；M2为积水区的像元均值；V为变异系数;Dw为相对差异。

1.3.2特征选取

变异系数反映了选取样本内像元统计特征值的离散程度，离散程度越小，即变异系数越小越容易根据该特征区分出该类地物。相对差异则反映的是地物之间在该统计特征上的差异程度，值越大表明越容易区分[19-20]。

本研究根据统计的变异系数和相对差异筛选出适于区分耕地和积水区的特征。首先根据耕地与积水区之间的27项相对差异值进行排序，选取相对差异大于70%的特征;然后再根据耕地的变异系数排序，选取变异系数小于30%的特征。最后选取色度(变异系数26.67%，相对差异73.33%)和绿色信息熵(变异系数20.59%，相对差异72.79%)作为耕地提取的最优特征，采用最大似然法进行提取，并统计耕地面积。

表2 图像特征统计

Table 2 Statistics of image features

指标耕地均值标准差变异系数/%积水区均值标准差变异系数/%相对差异/%饱和度145.2534.3523.65176.3991.1751.6921.44亮度0.390.1025.640.360.0925.007.69色度0.150.0426.670.040.0375.0073.33蓝色均值18.379.9554.1617.926.2334.772.45蓝色方差11.8714.96126.0311.2522.22197.515.22蓝色协同性0.490.2040.820.680.2435.29387.78蓝色对比度4.595.07110.464.237.69181.807.84蓝色相异性1.540.9762.991.061.21141.1531.17蓝色信息熵2.880.5820.142.021.0853.4729.86蓝色二阶矩0.080.0787.500.270.2488.89237.50蓝色相关性0.820.1821.960.600.3253.3326.83绿色均值21.836.2728.7216.735.5132.9323.36绿色方差6.217.74124.6410.2720.51199.7165.38绿色协同性0.550.1730.910.690.2434.7825.45绿色对比度2.742.85104.013.987.27182.6745.26绿色相异性1.170.6555.561.011.17115.8413.68绿色信息熵2.720.5620.590.741.11150.0072.79绿色二阶矩0.100.0770.000.300.2686.67200.00绿色相关性0.780.1823.070.580.3356.9025.64红色均值17.216.1735.8518.214.8326.525.81红色方差6.197.95128.438.9618.00200.9044.75红色协同性0.560.1730.360.690.2333.3423.21红色对比度2.722.88105.893.606.53181.3932.35红色相异性1.160.6656.900.971.09112.3716.38红色信息熵2.690.5821.561.941.0654.6427.88红色二阶矩0.100.0880.000.290.2482.76190.00红色相关性0.770.1924.680.530.3362.2631.17

1.4 无人机可见光植被指数

1.4.1植被指数选取

总结前人的研究，基于绿光通道与红光、蓝光通道间运算的以及植被指数之间运算的可见光植被指数见表3。

1.4.2植被指数计算与分析

分别计算各植被指数，得到植被指数空间分布如图2所示。总体上看，通过绿光通道与红光、蓝光通道间的运算，多数可见光植被指数增强了耕地信息的特征，EXG,VDVI,RGBVI及NGRDI的指数分布图中耕地呈现亮色，其它地物呈现暗色，而且色调越亮表示植被指数值越大;EXGR指数分布图中耕地与积水区均呈现亮色，说明EXGR指数并没有凸显耕地信息;WI指数分布图中耕地与道路、积水区颜色辨识度低，边界信息非常模糊，不利于耕地的提取。因此，根据可见光植被指数空间分布情况，选取空间分布色差明显的EXG,VDVI,RGBVI及NGRDI作为高潜水位煤矿区耕地信息提取的可见光植被指数。

表3 可见光植被指数

Table 3 Visible-band vegetation indexes

植被指数全称公式文献出处EXG过绿指数2G-R-B[21]EXGR超绿超红差分指数EXG-1.4R-G[22]WIWoebbecke指数(G-B)/(R-G)[23]VDVI可见光差异植被指数(2G-R-B)/(2G+R+B)[18]RGBVI红绿蓝植被指数(G2-RB)/(G2+RB)[24]NGRDI归一化绿红差异指数(G-R)/(G+R)[25]

注:R为红光波段反射率;G为绿光波段反射率;B为蓝光波段反射率。

阈值的确定是利用植被指数提取耕地信息的关键，即植被指数大于阈值的像素归为耕地。本文选取双峰直方图法确定各植被指数的阈值，经过植被指数灰度化后的图像，会在直方图上表现出2个明显的波峰，分别代表绿色植被和背景。选取2个波峰之间的波谷处作为阈值，可以很好地进行耕地的识别与提取。分别生成EXG,VDVI,RGBVI及NGRDI植被指数的灰度直方图，如图3所示，图中横坐标为各植被指数的值，纵坐标为像元个数，从图中可以看出，4种指数的直方图都具有明显的双峰特征，耕地提取的阈值分别为15.219 6，0.067 3，0.141 8，0.055 3。

图2 可见光植被指数分布

图3 植被指数直方图

2 结果与分析

2.1 不同方法提取结果分析

分别采用基于样本的面向对象、基于色彩与纹理特征、基于可见光植被指数3种方法开展高潜水位煤矿区沉陷耕地的提取，结果如图4所示。无人机影像中沉陷耕地呈现绿色，条带纹理明显。基于样本面向对象的沉陷耕地提取中，在多次试验的基础上，最终选择的分割与合并尺度分别为40，90，根据实地调查，选取耕地样本12个，采用支持向量机的方法完成沉陷耕地提取;通过色彩与纹理的图像增强，提取结果不仅反映了沉陷耕地的纹理特征，而且对于采煤沉陷导致的耕地稀疏和颜色的变化也有较为客观的表达;将4种植被指数的提取结果与研究区正射影像进行对比后发现，EXG指数的提取结果耕地颜色过分统一，并没有将部分耕地稀疏、颜色变化的特点准确表达，尤其是在沉陷积水区周围;NGRDI指数提取的结果中，积水区周围的部分耕地没有有效的提取出来。因此，EXG和NGRDI指数并不适用于高潜水位沉陷耕地的提取。而VDVI,RGBVI这2种可见光指数提取结果与实际情况较为接近，为了评价这2种指数的精度，在原始影像上建立耕地的感兴趣区域(ROI)，利用ROI与指数提取结果进行匹配，得出耕地总体分类精度，结果为VDVI指数(81.05%)高出RGBVI指数(71.38%)将近10个百分点，这说明VDVI是所选可见光植被指数中最适用于高潜水位煤矿区沉陷耕地的提取。

图4 沉陷耕地提取结果对比

2.2 耕地面积对比与分析

为了进一步比较两种沉陷耕地提取方法，利用验证区影像(图1(b))，以基于样本的面向对象的沉陷耕地提取结果作为参考值，统计了这2种方法提取的沉陷耕地面积，并得到了面积提取的误差，结果见表4。从表4可以看出，基于色彩与纹理特征的沉陷耕地提取面积与参考面积较为接近，误差较小。因此，仅就本研究而言，基于色彩与纹理特征法的精度高于VDVI可见光植被指数法。

表4 不同提取方法的结果与误差

Table 4 Results and errors of different extraction methods

参考面积/m2基于色彩与纹理特征面积/m2误差/%VDVI植被指数面积/m2误差/%2231.342079.166.81873.4116.0

3 讨 论

通过本研究可以看出，低成本、广泛使用的可见光相机获取的图像虽然是数码照片，缺少特定波段的光谱信息，譬如近红外波段，但是耕地在高空间分辨率遥感影像上呈现明显条带状纹理，通过对图像的进一步增强，基于色彩与纹理特征的提取方法要优于基于波段运算或者指数运算的可见光植被指数法，能够替代人工目视解译提取的耕地信息，解决了传统方法耗时长，投入工作量大，难以满足快速获取耕地信息的不足。虽然效果不如色彩与纹理特征提取法，但是可见光指数的构造方式在很大程度上增强了耕地信息，下一步对可见光植被指数的改进将是研究的重点。

作为煤矿区土地复垦或征地补偿工作中主要参考的边界，耕地损毁边界的准确划定这一关键问题至今没有得到很好的解决，主要原因是耕地损毁边界处的地面特征与非损毁区、非沉陷耕地相比无明显差别，边界具有模糊性。本研究不仅较为准确的提取了沉陷耕地信息，而且较为客观的反映了由于采煤沉陷导致耕地颜色、纹理、疏密等变化特征，这为高潜水位煤矿区沉陷耕地损毁边界的研究提供了基础。本研究使用的旋翼无人机起降方式简单，不需要跑道或弹射架，降落过程不需要降落伞，具有很高的安全性，在高潜水位煤矿区的数据采集中优势明显。但基于可见光波段的无人机遥感仍存在一定的局限性:① 电动无人机的续航时间较短，数据采集能力有限，飞行区域相对较小，更大区域数据的获取较为困难;② 无人机遥感是一种低空数据采集方式，获取的影像数据空间分辨率较高，可以达到厘米级，这对耕地农作物的精细识别提供了基础，但空间分辨率的提高，并不意味着影像自动解译精度的提高，BRUZZONE等[26]曾指出，遥感影像空间分辨率过高会导致类内方差变大，类间方差变小;空间分辨率过低导致混合像元增多，降低分类精度。因此，高分辨率无人机数据是否可以通过空间分辨率的降低来满足实际的需求还需要进一步研究;③ 可见光影像的波段数少，加之小型无人机缺乏定量光谱信息，一定程度上限制了地物的自动识别，同一地物在不同航片上可能存在不同的色彩灰度，这将导致计算得到的变异系数增大;④ 高分辨率遥感影像的信息提取多采用面向对象的图像技术，在沉陷耕地信息提取过程中，利用基于边缘监测的分割算法完成了影像的分割，但是分割尺度是否为最优尺度，还需要进一步的研究确定。

4 结 论

(1)通过比较高潜水位煤矿沉陷区典型地物对可见光波段的反射统计特性发现，耕地与积水区在3个波段均存在重合，需通过增强耕地的颜色及纹理特征，或者通过绿光通道与红光、蓝光通道间或植被指数间的运算，实现耕地的准确提取。

(2)对研究区无人机影像进行了色彩空间转换与二阶矩阵纹理滤波后，统计了耕地与积水区共27项色彩与纹理特征指标，根据均值和标准差计算了变异系数和相对差异值，最终选取色度和绿色信息熵作为耕地提取的最优特征，并采用最大似然法进行耕地提取。

(3)选取6种可见光植被指数作为沉陷耕地提取的备选指数，根据植被指数计算的结果分布图，选取了EXG，VDVI，RGBVI及NGRDI作为沉陷耕地提取指数，利用双峰阈值法确定了耕地提取阈值。比较提取结果得出，EXG和NGRDI指数无法全面、客观反映研究区实际情况，VDVI指数的耕地总体分类精度为81.05%，高于RGBVI指数的71.38%，是本研究中最适用于高潜水位煤矿区沉陷耕地提取的指数。

(4)以基于样本的面向对象的沉陷耕地提取面积作为参考值，通过比较验证区沉陷耕地提取面积及误差得出，基于色彩与纹理特征法的精度高于可见光植被指数法，更适用于高潜水位煤矿区沉陷耕地的提取。