基于无人机遥感图像提取农田微地形特征

杜蒙蒙，李瀚远，金 鑫，ROSHNIANFARD Ali，周 浩

（1.河南科技大学农业装备工程学院，河南 洛阳 471003；2.莫哈格达阿比德利大学农业与自然资源系，伊朗 阿尔达比勒 566199）

0 引言

近年来，极端暴雨天气与自然灾害频发，导致农田受损，影响农作物生长，降低了产量和品质[1-2]。在雨水冲刷或其他自然灾害等外界因素作用下，农田中往往出现尺寸不等的凸起与洼地等微地形特征，不利于进行农田耕作[3]。凸起特征在形态上表现为中间高、四周低的封闭地势，与经典地貌学中的山脊特征相似，指包含单个或多个高程极大值点的局部区域；而洼地特征则为四周高、中间低的封闭地势，与经典地貌学中的山谷相似，指包含一个或多个高程极小值点的局部区域[4]。

传统土地平整忽视了农田微地形特征的结构尺寸和位置分布，易造成平整地机具过载[5]。准确识别农田微地形特征，并在此基础上对其进行精准平整，可避免平整地作业过载率偏高的问题。数字地形模型（Digital Terrain Model，DTM）是分析地形特征的重要工具[6]。仲腾等[7]利用反转数字地形模型提取山顶特征，该方法需人为确定阈值，导致提取精度易受人员主观判断影响。易炜[8]采用图像处理算法划定最大标高值区域，根据人为设定的阈值提取山顶特征，该方法受人员主观判断的影响较大。DENG Y 等[9]使用多个地形因子综合评价用于提取山顶特征的阈值，在模型中有效剔除了不符合实际状况的山顶特征点。

本研究利用低空无人机遥感图像，建立高精度农田数字地形模型（Farmland Digital Terrain Model，FDTM），实现农田微地形特征的定位与轮廓提取，可为精准平整地作业提供方法支持，具有迫切的现实需求和重要的理论意义。

1 数据获取

DTM 空间分辨率大小直接决定了表达地形曲面的精细化程度，影响地形分析的精度。因此，建立高精度FDTM，可有效提高农田微地形特征的提取精度。研 究 区 域 位 于 陕 西 省 宝 鸡 市（107°9′14.893″E，34°39′182″N），其尺寸约为50 m×100 m。采用极飞C2000 型智能测绘无人机，搭载机械全局快门2 000 万像素的航拍相机，以距地面相对高度80 m 进行航拍测绘作业，共获取42 张航片，航片横向重叠率85%、纵向重叠率80%。基于SfM 技术，获取空间分辨率3 cm的FDTM，如图1 所示。试验田高程范围为853.527～854.821 m，高程极差为1.294 m。

图1 试验地块的FDTMFig.1 FDTM of experimental farmland

2 研究方法

2.1 确定局部窗口尺度

统计不同分析窗口尺度下FDTM 的高程方差，分析其随窗口尺度的变化规律，确定分析窗口的尺度范围，如图2 所示。

图2 高程方差变化率Fig.2 Change rate of elevational variance

由图2 可知，高程方差在[0，30] 像素和[60，100]像素区间内的变化率较大，呈发散状态；而在[30，60] 像素区间内的变化率较小，收敛于9.77×10−5m2。另外，相关文献表明，在高程方差收敛区间内，地形特征的识别精度随局部窗口尺度的增大而增加[10]。因此，本研究选择51 像素×51 像素窗口确定农田高程极大值点，共计315 个；结合31 像素×31像素、41 像素×41 像素和51 像素×51 像素3 个窗口，从农田高程极大值点集中定位凸起特征的顶点。

2.2 分析地形因子综合隶属度

农田现实环境中存在大量细碎土块和秸秆残茬，严重影响农田凸起特征的识别效率及精度。地形因子是数字地形分析的基础与核心内容，其主要分为宏观和微观地形因子[11]。本研究根据农田微地形结构，确定高程、坡度和地形起伏度作为农田高程极大值点的评价指标。高程表征某点沿铅垂线到大地水准面的距离；坡度计算如式（1）所示，指坡面竖直高度和水平距离的比值；地形起伏度计算如式（2）所示，指局部范围内地面高程最大值与最小值之差[11]。

式中S−坡度，（°）

fx−高程在水平方向的变化率

fy−高程在垂直方向的变化

RF−地形起伏度，m

Hmax、Hmin−局部窗口下的最大高程、最小高程，m

周侗等[12]研究表明，在单一窗口下受干扰因素影响，无法精准识别试验地块的微地形特征。因此，本研究首先利用式（3）对51 像素×51 像素窗口下的高程、地形起伏度和坡度数据进行归一化处理，使用熵权法对其进行权重赋值，如式（4）～式（6）所示，分别确定高程、地形起伏度和坡度的权重为0.384、0.487和0.129。于是，根据式（7）计算高程极值点在51 像素×51 像素窗口下的地形因子综合隶属度。进一步，分别计算31 像素×31 像素与41 像素×41 像素窗口下高程极值点的地形因子综合隶属度，并使用线性权重法确定31 像素×31 像素、41 像素×41 像素和51 像素×51 像素窗口的权重为0.252、0.333 和0.415，如式（8）所示。最终，利用式（9）计算多窗口地形因子综合隶属度，如表1 所示。

表1 多窗口地形因子综合隶属度Tab.1 Comprehensive membership degree of multi window terrain factors

式中i−地形因子序号，i=1，2，3，分别为高程、地形起伏度和坡度

j−高程极值点序号，j=1，2，…，315

k−窗口尺度大小，k=31，41，51

fijk−地形因子的归一化值

min（xijk）、max（xijk）−在窗口尺度k下第i个地形因子的最小值和最大值

xijk−在窗口尺度k下第j个高程极值点第i个

地形因子的值

Pijk−在窗口尺度k下第j个高程极值点第i个地形因子的概率

Hik−在窗口尺度k下各地形因子的信息熵

ωik−在窗口尺度k下各地形因子的权重

φjk−在窗口尺度k下3 个地形因子对各农田高程极值点的地形因子综合隶属度

γk−局部窗口的权重

ψjk−多窗口综合地形因子综合隶属度

3 结果与分析

3.1 提取农田凸起特征

随着多窗口地形因子综合隶属度的增大，相应的高程极大值点成为凸起特征顶点的概率随之增大。因此，多窗口地形因子综合隶属度的阈值取值过小，作物残茬或细小土块易被误判为农田凸起特征；反之，识别农田凸起特征的准确率降低。本研究利用斯特吉斯公式对多窗口地形因子综合隶属度进行抽样处理，以确定最佳的地形因子综合隶属度阈值[13]。通过式（10）和式（11）分别计算分级组数n=9 及分级组距D=0.075，确定多窗口地形因子综合隶属度与对应的高程极值点数量，绘制高程极值点数量的直方图，如图3 所示。

图3 多窗口地形因子综合隶属度频率Fig.3 Comprehensive membership degree frequency of multiwindow terrain factors

式中n−分级组数

C−农田高程极值点数量（315 个）

D−分级组距，m

φmax、φmin−多窗口地形因子综合隶属度的最大值和最小值

由图3 可知，在0.176～0.627 区间范围内，随着地形因子综合隶属度的增大，高程极值点的数量未见明显变化，表征农田中大量不规则分布的细碎土块与作物残茬；而在综合隶属度0.627 处高程极值点数量急剧减小，即出现尺寸显著增大的农田凸起特征。因此，确定0.627 作为地形因子综合隶属度的阈值，并提取表1 中地形因子综合隶属度大于该阈值的农田高程极值点，即农田凸起特征的顶点，共计16 个，实现农田凸起特征的定位，其空间分布如图4a 所示。

图4 提取农田凸起特征Fig.4 Extracting hump features of farmland

等高线无法定量描述坡度和地形起伏度等地形信息，但能够直观地表达地形结构特征。研究表明，等高距约等于地形图空间分辨率时可较为完整地表达地形信息[14]。因此，本研究利用3 cm 等高距的等高线图，并提取上述农田凸起特征顶点的等高线闭合区域，识别农田凸起特征的外形轮廓，如图4b 所示。

3.2 提取农田洼地特征

RDTM 是一种非客观存在的地形，以DTM 的平均高程作为水平面基准，经反向翻转构成的数字地形模型。RDTM 与DTM 存在反转对称关系，两者坡度相差180°，原始DTM 的高程值极小点表现为RDTM 的高程值极大点。因此，原始FDTM 的洼地特征在形态上表现为RFDTM 的凸起特征。建立研究区域的RFDTM，重复前文提取农田凸起特征的方法，获取9 个农田洼地特征，其空间分布如图5 所示。

图5 提取农田洼地特征Fig.5 Extracting concave features of farmland

4 结束语

以陕西省宝鸡市某试验区域为研究对象，采用极飞C2000 型智能测绘无人机进行航拍测绘作业，通过分析农田高程方差随局部窗口尺度的变化规律，确定适应于试验田微地形分析的局部窗口尺度范围为31～51 像素。利用相关系数确定高程、地形起伏度和坡度3 个地形因子，并计算农田高程极大值点在31 像素×31 像素、41 像素×41 像素和51 像素×51 像素窗口下的多窗口地形因子综合隶属度，根据阈值法定位16个凸起特征与9 个洼地特征，最终结合农田等高线图提取上述微地形特征的外形轮廓，实现农田微地形特征的定位与轮廓提取，为精准平整地作业提供支持。