0 引 言
煤炭资源作为我国的主要能源，在未来很长一段时间内其地位不会改变。 而煤炭开采势必会对土地、植被造成一定影响。 据估计，我国煤粮复合区的面积超过了总耕地面积的10%[1]。 特别是在高潜水位煤粮复合区，由于开采沉陷导致地下水位相对上升，部分农作物根系浸入水中，从而产量下降甚至绝产。 作物生长状况是确定土地损毁程度进而影响土地复垦难度和复垦方案选择的关键，也是矿山企业赔付农民的直接依据，因此，快速准确地获取采煤沉陷地作物绝产边界可更好地解决企地矛盾，也是实现高潜水位矿区精准治理的前提[2-7]。
传统的矿区沉陷地绝产边界都是靠实地勘察和测量来确定，费时、耗力，且由于缺少客观的标准易受人为因素影响。 使用无人机能快速获取影像信息，不仅提高了野外作业效率，且提供了直观的证据，可降低人为因素影响。 国内还未见使用无人机遥感影像来确定矿区沉陷地作物绝产边界的研究，不少学者在利用遥感影像提取地物边界方面做了相关研究[8-11]。 目前，无人机遥感[12]作为一种新型的遥感数据获取方式已经被广泛应用于各行各业[13]，如在农业中利用无人遥感来进行田间作物的分类[14]、作物长势的监测[15-16]、土地利用类型的变化[17]等，从而为精细化的农业管理提供服务。 利用遥感影像提取边界，关键是边缘信息的提取，即边缘检测，Canny 算子是进行边缘检测的经典算法[18]，但随着图像识别与边缘检测技术的不断深入发展，研究人员逐渐认识到在不同尺度分辨率下，边缘特征中包含的信息是不同的 [19-20]，因此，将最大类间方差法(大津算法，OTSU)与Canny 算子结合[21-23]，逐渐被广泛应用于图像分割[24-28]。 由于 OTSU 算法是对图像灰度值进行运算，且一般应用于离散数据，而作物绝产的主要特征是作物高程存在差异，且高程数据为连续值。 因此，不能直接采用OTSU 算法对其进行分割，需要进行改进。 基于此，笔者提出了一种改进的OTSU 算法，并将其用于无人机影像的绝产边界识别，最后以山东省济宁市某矿采煤沉陷区为例进行了方法验证。
1 研究区概况及数据源
1.1 研究区概况
研究区位于山东省济宁市境内某矿，地理坐标：116°50′49″E—116°56′56″E，35°24′11″N—35°31′25″N，是我国典型的高潜水位矿粮复合区，地面平均标高+50.72 m，开采煤层平均标高为-497.8 m，煤厚平均 8.30 m，倾角为 4°。 自 2014 年 8 月开采，地面沉陷最大达到7 m，已经形成了大面积的积水区，形成了典型的矿区积水盆地，沉陷区自内由外依次分布为芦苇、沼泽、滩涂、水淹的农作物。 该地区的主要种植农作物包括玉米和小麦，夏季降雨充沛，沉陷区的玉米更容易受到积水的胁迫，出现减产和绝产现象，故选择玉米为研究对象。 研究位置如图1 所示。
图1 研究区位置示意
Fig.1 Overview of research area
1.2 无人机数据获取及处理
无人机飞行时间为2017 年8 月12 日 11：00—11：30，天气晴朗，室外温度 33 ℃左右，南风 2 级。采用四旋翼大疆经纬100 无人机搭载Zenmuse X3数码相机，设定飞行高度为128 m，航速为9 m/s，航向重叠度和旁向重叠度分别为80%和60%。
综合考虑像控点的布设原则以及研究区的实际情况，布设8 个像控点(图1)，并用RTK 测定三维坐标。 利用Pix4Dmapper 软件对无人机获取的影像进行数据处理，通过特征点匹配、像控点校正、相片拼接和空三加密后形成研究区的正射影像图(DOM)和数字表面模型(DSM)，经计算空间分辨率约为5.5 cm。
2 OTSU算法原理及其改进
2.1 最大类间方差(OTSU)阈值分割算法
最大类间方差(OTSU)阈值分割算法是一种使类间方差最大的自动确定阈值的方法，该方法具有简单、处理速度快的特点，是一种常用阈值分割方法。 基本思想如下：设图像像素的总数为N，灰度区间为[0，L - 1]，其中L 为灰度级总数，对应灰度级i 的像素数为 ni(i ＝ 0，1，2，…，L -1) ，灰度级 i 出现的概率pi 为
将图像中像素按灰度使用阈值T 分成2 类C0和 C1， C0 是由灰度在 [0，T]之间的像素组成， C1 是由灰度在[T + 1，L - 1]之间的像素组成。 对于灰度分布概率，整幅图像的均值为
而类C0 和C1 的均值分别为
其中：w0 和w1 分别为被分到C0 类和C1 类的概率。
由式(3)—式(5)可得
于是，类间方差为
利用式(5)和式(6)，由式(7)得
再利用式(6)和式(8)得
让T 在[0，L - 1]范围内依次取值，使类间方差
最大的T 值即OTSU 算法的最佳阈值。
2.2 OTSU阈值分割算法的改进
由以上OTSU 算法原理分析可知，采用传统的OTSU 算法只能针对离散的像素灰度来确定最佳的分割阈值，而本文研究对象DSM 的像素为非离散的，因此需要予以改进。 改进算法的核心包括2 个方面：①由于研究区DSM 影像的每个像素单元记录的是高程，是连续的，而OTSU 一般只适用于离散变量的情况，所以不能使用OTSU 算法直接对研究区DSM 影像进行分割，需要将高程进行离散化改进。②由图2 可以看出，从研究区DSM 影像中的沉陷区由内而外任取一条剖面线，根据无人机点云数据提取作物高点高程，可见，在绝产边界内、外作物高程上有很大的差异，因此在提取绝产区和非绝产区时，可以利用作物高程上的差异对DSM 影像进行分割，即用高程代替OTSU 算法中的灰度，进行图像分割。
图2 作物绝产边界两侧的作物高程差异
Fig.2 Elevation difference on both sides of crop yield boundary
改进的OTSU 算法具体操作如下：
1)高程离散化：假设DSM 影像的最小高程和最大高程分别为Dmin和Dmax，选取一个初始的高程分割份数N0，则对应的初始高程间隔S 为
以S 为间隔在Dmin和Dmax之间取值。
2)利用原有的OTSU 算法，寻找N0 下的最大类间方差σ0 和对应的高程h0，同时选新的高程分割份数 Ni ＝ 2N0。
3)计算新的高程分割份数Ni 下的高程间隔Si：
重复步骤2 的过程，得到新的类间方差σi 及其对应的高程hi 。
4)当hi - hi-1 小于给定的阈值时，则停止迭代。
3 农作物绝产边界检测与验证
3.1 基于改进OTSU算法的DSM分割
经过改进后的OTSU 可适用于本研究，将数据输入迭代完成即可求出高程阈值，按该阈值将图像二值化分割，图像分割及二值化过程如下。 从研究区数字表面模型(DSM，图 3a)获取 Dmin ＝45.33，Dmax ＝48.27，选择初始分割份数 N0 ＝ 10，2 次运算的高程之差为hi - hi-1 ，将hi - hi-1 阈值设置为0.01，利用改进的OTSU 算法，经过4 次迭代，得到高程分割阈值hi ＝47.40 m。 然后，将高程大于47.4 m 的区域确定为非绝产作物区域，像素置为白色，反之确定为绝产区域，置为黑色，得到二值图像(图3b)。
3.2 作物绝产边界的检测
边缘检测的基本思想是，利用图像在边缘像素点周围的灰度存在明显的变化，通过边缘检测算子量化该变化的大小，从而达到提取边缘的目的。 在上述图像分割使灰度图像尽可能地简化，从而使得图像中的信息更加纯粹，边缘值变化地更加明显，有利于边缘检测的准确性。 为了得到二值图像中黑白相交边界信息，使用Canny 算子进行边缘检测，提取到的绝产边界如图3c 所示。
图3 图像分割及边缘检测
Fig.3 Image segmentation and edge detection
3.3 绝产边界的验证
试验中无人机影像可达厘米级别的空间分辨率，在研究区的DOM 影像图上作物绝产边界表现出就是一条由正常作物生长区域向不正常作物生长区域的过渡，人眼较容易区分，如正常生长玉米与非正常生长玉米、杂草等在影像上色调差异明显，因此采用目视法进行验证评价。 在研究区正射影像图(DOM，图4)上通过目视选取12 个验证样本，其中11 个样本中包含边界线，8 个正确。 图4 为其中4个典型的验证样本位置及边界提取情况，图中红线表示提取的边界。 包括玉米明显死亡点(图4c)、玉米生长状况受影响点(图4a)以及玉米已经死亡并生长出杂草等(图4b、图4d)。 由验证样本可以看出，正常生长玉米与枯萎、受影响生长不良、杂草等在高程上存在明显差异，通过高程可以较为清晰地区分，而改进的算法可以用于确定高程阈值的图像分割，提取作物绝产边界的结果比较准确。
图4 农作物绝产边界的验证样本
Fig.4 Validation samples of crop yield boundaries
4 结 论
在高潜水位煤粮复合区，采煤沉陷地由于地下水位的相对上升，会造成耕地作物的根部被水浸泡，从而导致作物绝产。 快速准确地获取作物绝产边界是高潜水位煤粮复合区环境综合治理和沉陷地复垦工作的前提。 利用采煤沉陷地的无人机遥感影像中绝产边界两侧作物在高程上存在较大差异的现实，改进了OTSU 图像分割算法，提出了一种基于影像的作物绝产边界识别方法。 主要研究结论如下：
1)提出了改进的OTSU 图像分割算法，实现了采煤沉陷区数字表面模型(DSM)中高程数据的离散化，使其能够应用于OTSU 图像分割算法。
2)改进的OTSU 图像分割算法科学利用了作物绝产边界两侧的作物高程差异，拓展了OTSU 图像分割算法的应用。
3)通过研究区正射影像图(DOM)的验证样本分析，该方法识别的作物绝产边界较为准确。
4)研究的缺陷表现在识别研究区作物绝产边界时仅从研究区的DSM 上入手，没有结合DOM 影像上的颜色特征，若两者结合使用可能会达到减少部分绝产边界的偏移；此外，本研究没有考虑到结合光谱植被指数如NDVI 去反映研究区作物生长状况，这也是后续需要研究的重点。
[1]李 晶，胡振琪，李立平.中国典型市域煤粮复合区耕地损毁及其影响[J].辽宁工程技术大学学报：自然科学版，2008(1)：148-151.LI Jing，HU Zhenqi，LI Liping.Farmland damage and its impacts in overlapped areas of cropland and coal resources in one typical city of China[J].Journal of Liaoning Technical University：Natural Science，2008(1)：148-151.
[2]胡振琪，赵艳玲，程玲玲.中国土地复垦目标与内涵扩展[J].中国土地科学，2004(3)：3-8.HU Zhenqi， ZHAO Yanling ，CHENG lingling.Extension of goal and meaning of land reclamation in China [J].China Land Science，2004(3)：3-8.
[3]卞正富，雷少刚，金 丹，等.矿区土地修复的几个基本问题[J].煤炭学报，2018，43(1)：190-197.BIAN Zhengfu， LEI Shaogang， JIN Dan， et al.Several basic scientific issues related to mined land remediation[J].Journal of China Coal Society，2018，43(1)：190-197.
[4]桑李红，付梅臣，冯洋欢.煤矿区土地复垦规划设计研究进展及展望[J].煤炭科学技术，2018，46(2)：243-249.SANG Lihong， FU Meichen， FENG Yanghuan.Progress and prospect of research on land reclamation planning and design in mining area[J].Coal Science and Technology，2018，46(2)：243-249.
[5]胡振琪.我国土地复垦与生态修复30 年：回顾、反思与展望[J].煤炭科学技术，2019，47(1)：25-35.HU Zhenqi.The 30 years’land reclamation and ecological restoration in China：review，rethinking and prospect[J].Coal Science and Technology，2019，47(1)：25-35.
[6]胡振琪，魏忠义，秦 萍.矿山复垦土壤重构的概念与方法[J].土壤，2005(1)：8-12.HU Zhenqi，WEI Zhongyi，QIN Ping.Concept of and methods for soil reconstruction in mined land reclamation[J].Soils，2005(1)：8-12.
[7]李 晶，刘喜韬，胡振琪，等.高潜水位平原采煤沉陷区耕地损毁程度评价[J].农业工程学报，2014，30(10)：209-216.LI Jing，LIU Xitao，HU Zhenqi，et al.Evaluation on farmland damage by underground coal-mining in plain area with high groundwater level[J].Transactions of the Chinese Society of Agricultural Engineering，2014，30(10)：209-216.
[8]齐文娟，杨晓梅.江西省都昌县山地与平原植被界线提取[J].地球信息科学学报，2017，19(4)：559-569.QI Wenjuan， YANG Xiaomei.Mountain and plain vegetation boundaries extraction in Duchang County Province Jiangxi[J].Journal of Geo-information Science，2017，19(4)：559-569.
[9]李 晶，焦利鹏，申莹莹，刘乾龙.基于IFZ 与NDVI 的矿区土地利用/覆盖变化研究[J].煤炭学报，2016，41(11)：2822-2829.LI Jing，JIAO Lipeng，SHEN Yingying，et al.Land use and cover change in coal mining area by IFZ and NDVI[J].Journal of China Coal Society，2016，41(11)：2822-2829.
[10]怀保娟，李忠勤，孙美平，等.多种遥感分类方法提取冰川边界探讨：以喀纳斯河源地区为例[J].干旱区研究，2013，30(2)：372-377.HUAI Baojuan，LI Zhongqin，SUN Meiping，et al.Discussion on RS methods for glacier outline detection：a case study in headwaters of the Kanas Rive[J].Arid Zone Research，2013，30(2)：372-377.
[11]李亚平，苗 放，黄舒寒懋.利用塌陷积水对采煤塌陷区进行监测：以开滦煤矿为例[J].内蒙古石油化工，2009，35(5)：1-3.LI Yaping，MIAO Fang，HUANG Shuhanmao.Kailuan coal mine collapse，as well as the causes of stagnant water in time and space to monitor the role of research[J].Inner Mongolia Petrochemical Industry，2009，35(5)：1-3.
[12]李德仁，李 明.无人机遥感系统的研究进展与应用前景[J].武汉大学学报：信息科学版，2014，39(5)：505-513，540.LI Deren，LI Ming.Research advance and application prospect of unmanned aerial vehicle remote sensing system[J].Geomatics and Information Science of Wuhan University，2014，39(5)：505-513，540.
[13]金 伟，葛宏立，杜华强，等.无人机遥感发展与应用概况[J].遥感信息，2009(1)：88-92.JIN Wei，GE Hongli，DU Huaqiang，et al.A review on unmanned aerial vehicle remote sensing and its application[J].Remote Sensing Information，2009(1)：88-92.
[14]刘建刚，赵春江，杨贵军，等.无人机遥感解析田间作物表型信息研究进展[J].农业工程学报，2016，32(24)：98-106.LIU Jiangang，ZHAO Chunjiang，YANG Guijun，et al.Review of field-based phenotyping by unmanned aerial vehicle remote sensing platform[J].Transactions of the Chinese Society of Agricultural Engineering，2016，32(24)：98-106.
[15]陈 雯，武 威，孙成明，等.无人机遥感在作物监测中的应用与展望[J].上海农业学报，2016，32(2)：138-143.CHEN Wen， WU Wei， SUN Chengming， et al.Application and prospect of UAV remote sensing in crop monitoring[J].Acta Agriculturae Shanghai，2016，32(2)：138-143.
[16]肖 武，陈佳乐，赵艳玲，等.利用无人机遥感反演高潜水位矿区沉陷地玉米叶绿素含量[J].煤炭学报，2019，44(1)：295-306.XIAO Wu， CHEN Jiale， ZHAO Yanling， et al.Identify maize chlorophyll impacted by coal mining subsidence in high groundwater table area based on UAV remote sensing[J].Journal of China Coal Society，2019，44(1)：295-306.
[17]韩文霆，郭聪聪，张立元，等.基于无人机遥感的灌区土地利用与覆被分类方法[J].农业机械学报，2016，47(11)：270-277.HAN Wenting，GUO Congcong，ZHANG Liyuan，et al.Classification method of land cover and irrigated farm land use based on UAV remote sensing in irrigation [J].Transactions of The Chinese Society of Agricultural Machinery，2016，47(11)：270-277.
[18]CANNY J.A computational approach to edge detection[J].IEEE T-PAM I，1986，8(6)：679-698.
[19]MALLAT Stephane.Zero-crossings of wavelet transform [J].IEEE Tans.Information Theory，1991，37(4)：19-1033.
[20]MALLAT Stephane，WEN Liang Hwang.Singularity Detection and Processing with Wavelets [J].IEEE Trans.On Information Theory，1992，38(20)：617-643.
[21]王 植，贺赛先.一种基于Canny 理论的自适应边缘检测方法[J].中国图象图形学报，2004(8)：65-70.WANG Zhi，HE Saixian.An adaptive edge-detection method based on canny algorithm [J].Journal of Image and Graphics，2004(8)：65-70.
[22]何新英，王家忠，孙晨霞，等.基于数学形态学和Canny 算子的边缘提取方法[J].计算机应用，2008 (2)：477-478，483.HE Jiaying，WANG Jiazhong，SUN Chenxia，et al.Edge detection method based on mathematical morphology and canny algorithm[J].Computer Applications，2008(2)：477-478，483.
[23]李华强，喻擎苍，方 玫.Canny 算子中Otsu 阈值分割法的运用[J].计算机工程与设计，2008(9)：2297-2299.LI Huaqiang， YU Qingcang， FANG Mei.Application of Otsu thresholding method on Canny operator[J].Computer Engineering and Design，2008(9)：2297-2299.
[24]王 楚，王亚刚.基于拟合的亚像素边缘检测的算法比较[J].软件导刊，2018，17(2)：60-62，67.WANG Chu，WANG Yagang.Comparison of algorithms for subpixel edge detection based on fitting[J].Software Guide，2008，17(2)：60-62，67.
[25]张 超，王志浩，杨建宇，等.基于Canny 算子的农田线状工程地物自动提取方法[J].农业机械学报，2015，46(2)：270-275.ZHANG Chao，WANG Zhihao，YANG Jianyu，et al.Farmland linear project feature auto-extraction method based on Canny algorithm[J].Transactions of The Chinese Society of Agricultural Machinery，2015，46(2)：270-275.
[26]景晓军，李剑锋，刘郁林.一种基于三维最大类间方差的图像分割算法[J].电子学报，2003，31(9)：1281-1285.JING Xiaojun，LI Jianfeng，LIU Yulin.Image segmentation based on 3-D maximum between-cluster variance[J].Acta Electronica Sinica，2003，31(9)：1281-1285.
[27]NG HF.Automatic thresholding for defect detection[J].Pattern Recognition Letters，2006，27(14)：1644-1649.
[28]WEN J，LI S，SUN J.A new binarization method for non-uniform illuminated document images[J].Pattern Recognition，2013，46(6)：1670-1690.
