干旱区水浇地提取研究

2021-07-01顾云志曹月娥赵珮珮

江苏农业科学 2021年10期

顾云志曹月娥赵珮珮

摘要：以新疆和田地区墨玉县阔依其乡萨达克村为研究区，利用高分二号遥感影像的纹理和光谱特征，采用阈值分割与数学形态学算法，研究出一种面向对象的多特征融合的干旱区水浇地地块提取方法。结果表明，利用阈值分割算法和数学形态学算法4个几何特征参数面积、周长、完整度、延伸率可以有效剔除细小图斑，准确提取水浇地边界及面积。通过研究区实地核查，地块面积提取精度为90.08%，利用随机样本点验证水浇地总体分类精度达到96.33%，Kappa系数为91.72%。

关键词：水浇地;地块提取;特征融合

中图分类号： S127 文献标志码： A 文章編号：1002-1302（2021）10-0176-07

水浇地作为中国北方地区粮食生产的主要载体，掌握其使用现状及面积变化，对于保证国家粮食安全和农业可持续发展具有重要作用。以往水浇地监测是依据地面调查方法，进行逐级汇总，对于小范围地域的精细化种植管理具有重要作用[1]，但这种方法消耗大量的人力、物力及财力，实施效率低下。随着农业现代化的发展，大范围实时耕地监测成为一种新趋势，航空摄影、遥感卫星提供高时空分辨率图像已经成为耕地监测的主要数据源[2-3]。利用图像通过地物的光谱、纹理、形状特征能够比较准确地提取水浇地的分布现状，能够对农作物的播种、生长状况进行监测，同时对粮食估产也可以提供基础数据。

当前，多光谱、高光谱和光学合成孔径雷达（SAR）等图像数据已经在资源调查、环境监测等方面广泛应用[4]，多光谱和高光谱数据融合成为提取土地利用的主要手段。根据研究内容及尺度，现有研究数据源涉及低、中、高等3种分辨率图像，其中MODIS[5]、NOAA等低空间分辨率的影像主要用于大尺度的植被变化监测及生态环境变化研究，而中高空间分辨率影像Sentinel-2、LandSet TM、SPOT、QuickBird、WorldView-1等在中小尺度地物识别及土地利用类型变化监测中应用较多。随着中国高分遥感系列卫星取得重大成功，其高时空分辨率影像产品应用于建筑物、河流湖泊、农业种植等各类地物提取[6]，为土地合理利用开发提供了基础数据保障。水浇地提取关键就是对图像进行分割，主要利用图像灰度值特征[7]、纹理特征及形态学特征[8]对地块进行分割。图像分割算法常用的主要有4个：基于阈值的分割方法[9-10]，利用二值化图像灰度直方图进行分割[11];基于边缘检测的分割方法[12]，通过计算图像局部的微分算子进行分割[13]，主要算法包括Roberts算子、Sobel算子、Canny边缘检测器[14];基于区域的分割方法，依据像素的灰度级进行分割与合并[15];形态学分水岭算法，首先将图像抽象化为三维，再依据图像梯度划定其分水岭分割线[16]，进而再对图像分割。基于同一图像高频和低频分量的不同，在不同地物边缘分割时采用双阈值和多结构数学形态学算法进行[17]。深度学习蓬勃发展过程中，聚类分析、模糊集理论[18]、基因编码、小波变换[19]等理论逐步融入到图像分割中，使分割精度提升。如何利用图像自身特征及适宜算法提升耕地提取精度成为农业信息获取的研究热点。

新疆南疆地区土地利用状况复杂，防风沙林地及土地扩张，导致耕地破碎化程度变高。因此，为提升耕地利用效率需要准确提取耕地的边界和面积。本研究利用GF-2遥感影像，选取新疆和田地区墨玉县萨达克村为研究区，以萨达克村耕地（水浇地）为研究对象，综合利用影像光谱、纹理、形状特征，结合阈值分割、数学形态学算法，研究一种面向对象的多特征融合的水浇地提取方法，以期为新时期干旱区农业发展提供技术支持。

1 材料与方法

1.1 研究区概况

萨达克村位于新疆和田地区墨玉县（图1），处于喀拉喀什河冲积扇平原、塔克拉玛干沙漠南缘，属于暖温带干燥荒漠气候，常年受风沙侵害。萨达克村经济发展以第一产业为主，粮食作物耕种是耕地利用的主要形式，其土地利用现状为耕地与居民点、交通用地、防护林交错，造成耕地信息统计困难。

1.2 数据来源及预处理

本研究采用2019年12月5日国产高分二号（GF-2）影像数据，GF-2卫星拥有极高的空间分辨率，对于地物类型提取拥有高精确度，其有效荷载技术指标见表1。

基于研究目的及数据特点，利用几何校正模型（RPC）对影像数据进行正射校正，然后依据中国资源卫星应用中心公布的高分二号的定标参数和光谱响应函数对处理后的影像进行辐射校正，并采用基于最小二乘法的光谱锐化（Gram-Schmidt pan sharpening，GS）影像融合方法对全色影像与多光谱影像进行叠加融合，同时对图像进行去噪处理，并利用行政边界进行裁剪，获取萨达克村的GF-2影像数据。

1.3 研究方法

1.3.1 纹理特征提取纹理信息是地表覆被类型细节的反映，不同类型地物的纹理特征表现形式是不同的[20]，主要通过图像区域内灰度变化为特征，显示为色调作为等级函数在空间上的变化。图像纹理特征通过交互式直方图均衡化拉伸实现辐射增强，通过凸显地物的内部具体纹理构造，对不同地物进行细致划分。由于影像各波段纹理特征具有同一性，因此利用单波段对影像纹理特征提取就能够对地物结构、影纹等纹理基本单元进行深度数据挖掘[21]。图像多光谱、纹理特征两者结合可以体现图像中不同地物类型的宏观属性与微观特征，有效提升整个图像中的地物识别度及分类精度。因此，本研究利用红波段对水浇地纹理特征进行提取，采用灰色空间相关性矩阵（gray level co-occurrence matrix，GLCM）计算图像纹理特征值。其中，纹理特征值是通过计算图像中特定距离和方向计算某个像元周边各个灰度级之间的联合概率密度获取该像元像素的最大概率灰度值，常用的主要纹理特征因子[22]有均值ME=∑N-1i，j=0iPi，j、方差VA=∑N-1i，j=0iPi，j（i-ME）2、对比度CO=∑N-1i，j=0iPi，j（i-j）2、同质性HO=∑N-1i，j=0iPi，j1+（i-j）2、异质性DI=∑N-1i，j=0iPi，j|i-j|、熵EN=∑N-1i，j=0iPi，j（-lnPi，j）、二阶矩SM=∑N-1i，j=0iPi，j2、相关性CC=∑N-1i，j=0iPi，j（i-ME）（j-ME）VAiVAj。其中Pi，j=Vi，j/∑N-1i，j=0Vi，j，式中，Vi，j为窗口单元格的第i行和第j列的数值;N为行列数。

纹理信息作为空间范畴内特征，选择合适的观察窗口对地物边界提取十分重要。较小的观察窗口利用波段纹理同一性的变化可以区分不同地物类型纹理边界，但图斑细碎化程度高;较大的观察窗口在纹理同一性检测方面效果好，能够准确识别连片地物，但纹理边界辨别度较差，在不同地物分界线划定时误差率高。因此本研究采取多个像元观察窗口（3×3、5×5、7×7、9×9）提取图像纹理特征，然后将不同观察窗口纹理提取结果进行对比，确定最适宜观察窗口，以准确分离出各地类边界。

1.3.2 阈值分割及边缘检测高分二号遥感影像中，同一地物中内部光谱纹理信息相对一致，但其中普遍存在“椒盐”现象，呈现出由噪声产生的各种形状不一、破碎的像元，严重影响提取准确性。通过阈值分割可以有效消除噪声影响，但在提取过程中，水浇地与小块建设用地、林地图斑相似性较高，会出现较多相互融合像元，不易被分割，因此采用Canny边缘检测算子对地物边界进行分割提取。Canny边缘检测算子[23]主要分为5个步骤：利用高斯滤波器将图像平滑及消除噪声;计算像元区域周围梯度强度和近似偏移方向;通过非极大值抑制消除像元杂散响应;利用双阈值或多阈值探测潜在的地类边缘;抑制孤立的像元弱边缘，最终完成整个图像各地类的边缘检测。Canny边缘检测算子能够从不同视觉方向提取地物类别结构信息，准确分离每个地类的边界，将影像细分为不同地物图斑，对于干旱区水浇地中破碎地块提取效果显著。

1.3.3 数学形态学算法地理环境因素条件下，干旱区水浇地四周环绕着大片林地，因此在阈值分割和边缘检测之后还需要利用数学形态学算法对地类与已定义的目标结构元素进行特定运算，以得到水浇地的具体形态边界，实现地块识别与分类[24]。形态学是基于特定目标形状的一种算法，体现为图像中每个像元的分类结果值取决于邻近区域的像元像素值及整体目标的形态。数学形态学运算中，主要是利用结构采集获取图像的信息，通过自身所拥有的灰度、形状特征，在移动过程中对其所经过的区域进行逻辑运算，并将运算的结果保留在原始的像元点上[25]。目标结构元素中形状特征的类型主要有任意型、线型、方型、矩型、菱型、碟型等。对于干旱区水浇地而言，其呈现为不规则形状，存在显著的向异性特点，因此选用任意型结构元素。任意型结构能够有效地去除不同尺度下的噪声和与目标不相关的图像像元，最终留下水浇地的图斑信息。

数学形态学基本运算包括膨胀、腐蚀、开运算、闭运算等[26]。膨胀运算能够将目标地类背景融合，填充目标中所存在的孔隙;腐蚀运算使目标缩小，可以消除图像中小于结构元素的像元;开运算在膨胀和腐蚀运算的基础上，突出图像明亮细节及平滑图像的较大地类的边缘;闭运算基于膨胀及腐蚀运算突出图像中较暗的细节，并使地类的位置及面积保持不变。基于数学形态学运算类型特征，为准确获取水浇地的边界及面积信息，本研究采取开运算与闭运算相结合的混合算法提取水浇地信息。

1.3.4 混淆矩阵精度验证水浇地提取精度结果利用图像分类常用的混淆矩阵进行验证。通过将实地调查结果与提取结果进行混淆矩阵运算，论证多特征融合算法的准确性与可靠性。混淆矩阵的精度验证一般考虑总体分类精度和Kappa系数（K）[27]，两者数值越高，则表明采用的分类方法越好，分类效果越显著。

总体分类精度=（TP+TN）/（TP+FN+FP+TN）;K=N∑xk-∑xkk∑x∑kN2-∑xkk∑x∑k

式中：TP为正确分类到水浇地的像元数，为真正类;TN为准确分到非水浇地的像元数，为真负类;FN为水浇地误分到其他类别中的像元数，为假负类;FP为其他类误分到水浇地像元数，为假正类。N=TP+FN+FP+TN，代表像元总数;x为类别数目，k为对角线上的数目，∑xkk为混淆矩阵对角线上的和，∑xkk∑x∑k为某一类别中真实参考像元总数与该类别中被分类的像元总数之积对所有类别求和的结果。

2 结果与分析

2.1 影像纹理特征提取最优窗口选择

本研究基于ENVI中的Co-occurrence 度量工具，采用3×3、5×5、7×7、9×9等4种不同观察窗口确定最优窗口。利用最优观察窗口初步划分水浇地与其他地类的边界，同时弱化其他用地对水浇地提取的影响。通过卷积滤波中高通滤波器消除特定空间频率，将图像做增强处理，使其能够保持图像的高频信息，锐化地块纹理、边缘信息。通过基于GLCM的主要纹理特征因子组合得出纹理信息，对不同观察窗口进行分析。各个观察窗口纹理特征顯示及提取结果见图2。结果表明，3×3像元的观察窗口较小，会使图斑细碎化程度变高，水浇地边缘区域像元被误分到道路、房屋等建设用地中，同时还造成水浇地边缘过于曲折;随着观察窗口变大，7×7、9×9像元观察窗口中图像纹理特征值计算量变大和计算时间延长，使得较小的农业设施及道路被同化到水浇地中，体现为水浇地边界向外延伸，突破固有水浇地的形状，提取效果变差。从整个图像纹理信息提取效果来看，随着像元窗口增大，光谱反射率高的建筑用地信息更加完整，而反射率较低的林地、水浇地信息则比较破碎。通过分析不同地类内部纹理特征，分布散乱的林地信息及形状规则的水浇地信息被从整个图像中分离出来，形成单独的纹理图层。经过平衡地类误分割概率、计算时间及计算量，最终选取5×5观察窗口作为最优观察窗口提取地物信息。

2.2 纹理与光谱特征结合

利用图像红波段的光谱特征提取水浇地信息，分类结果显示其与建设用地分离程度明显，但水浇地与林地仍存在较多像元混合的现象（图3）。从整个过程来看，纹理信息特征展现了地类内部细节特点，而光谱特征依据每个地类的光谱响应曲线能够判定其具体属性，将图像2个特征叠加后提取的水浇地结果见图4。结果显示，单波段提取水浇地过程中，水浇地与成片建设用地可以分离，但会出现与小面积建筑物、道路相互混淆的情况，而且地块完整度较低，误分概率大。但通过纹理与光谱信息结合则可以剔除细碎建筑用地及分散林地的干扰，提高水浇地地块完整度，更好地凸显水浇地边缘信息。

2.3 阈值分割与形态学算法提取水浇地

基于纹理及光谱特征信息提取的结果，通过ENVI中基于特征提取的分割（segment only feature extraction workflow）工具对分割阈值进行优化调整，基于IDL平台利用Canny边缘检测算法对水浇地边界信息进行提取（图5），主要代码如下。结果显示，建设用地周边的破碎地块与林地仍存在部分像元相连，需要进一步分割。

;Read a greyscale image

image=READ_TIFF（FILEPATH（‘Area_clip55_canny.tif，$ SUBDIRECTORY=[‘examples，‘data]））

;Resize the image

image=REBIN（image，505，605）

image=reverse（image，2）;数组x;水平翻转：reverse（x，1）;垂直翻转：reverse（x，2）

;Perform edge detection using defaults

filtered=CANNY（image）

;Create a window

WINDOW，XSIZE=758，YSIZE=605;

;Display the original and filtered images

TVSCL，image，0

TVSCL，filtered，1

对于水浇地中破碎地块与分散林地像元混合的现象，采用数学形态学算法进行细分割。通过对阈值分割后的影像进行数学形态学开运算和闭运算，锐化水浇地边缘信息和建设用地周边地物的细节信息，再利用水浇地固有的几何特征（面积Area、周长Perimeter、完整度Solidity、延伸率Elongation）去除非水浇地图斑。经过多次几何特征参数调整，使水浇地提取达到良好效果。本研究采用的几何特征参数为：Ar<40，去除孤立的细碎图斑;Ar>40、Per<30去除规则的细小图斑;Sol>0.25、El<1去除其他不規则图斑。通过上述算法及几何特征运算，水浇地最终提取结果见图6。

2.4 精度验证

通过ArcGIS将水浇地提取结果转换成矢量数据，并对图层进行投影坐标转换，利用获取的水浇地分布状况，统计水浇地面积。采用混淆矩阵中总体分类精度及Kappa系数对分类精度进行评价。为保证精度准确性，通过在区域内随机生成300个样本点（图7），获取样本点落在水浇地（图8）及其他用地类型上的数量，进而得到TP、TN、FP、FN数据，制作混淆矩阵（表2）。基于多特征融合算法提取结果与实地调查结果汇总比较，得到最终分析结果。精度评价结果显示，面向对象的多特征融合算法提取水浇地面积为97.29 hm2，而实地核查统计面积为108.00 hm2，自动提取精度达到90.08%，能够准确提取出水浇地空间位置分布及其现有面积。通过随机生成的样本点评价分类精度，结果显示水浇地总体分类精度达到96.33%，Kappa系数为91.72%，证明面向对象的多特征融合算法能够有效满足水浇地提取需求，能够为干旱区农业可持续发展提供基础耕地信息支持。

3 结论与讨论

本研究利用GF-2影像的纹理、光谱、几何特征，结合阈值分割、数学形态学算法，研究一种面向对象的多特征融合的干旱区水浇地提取方法。通过验证该算法能够准确地提取干旱区水浇地边界及面积，经过与实地核查水浇地面积数据对比，提取精度达到90.08%。通过随机样本点验证，结果显示水浇地总体分类精度达到96.33%，Kappa系数为91.72%，地块提取精度较高。在以往研究中，对于耕地提取主要是利用影像时相特征并结合光谱、纹理特征进行[28-29];研究区域主要集中在农业发展条件良好、耕地平整地区;提取方法主要集中在图像物理分割的算法方面[30]。而本研究依据已有研究成果，将几个方面有机融合，提出一种基于多特征融合的水浇地提取方法，为准确提取水浇地边界及面积作出参考。此提取方法能够有效发现耕地的变更状况，准确了解水浇地的分布状况，为新时期干旱区农业发展提供基础耕地利用数据。

研究发现，水浇地与林地混交在道路及居民点附近，是造成水浇地提取困难、影响提取精度的主要因素。这种现象广泛分布于中西部干旱区，主要是由于以耐旱树种为代表的林地具有重要的生态作用，林地为保障以耕地及建设用地为代表的生产生活用地正常运转被种植在其周边区域。农林混交种植一直被看作农业信息提取精度提高的最大阻碍，随着高分辨率影像作为数据源，并结合高程、蒸散发等多源数据和差值、土壤调节、比值等植被指数[22]，耕地信息提取精度才得到极大提高。未来耕地信息获取必将以多源数据为基础，依据长时间序列的图像光谱、纹理、几何特征，并结合多种算法及模型进行。

本研究基于图像纹理特征和光谱特征，运用阈值分割及数学形态学算法，并结合水浇地几何特征的面积、周长、完整度、延伸率等4个参数对图像细分割，最终自动提取获得水浇地边界及面积。在研究过程中，阈值分割及几何特征中参数设置都是基于过往经验及经过人为因素多次调试获得，以寻求最佳参数值，获取最优提取结果。最优参数值获取需要在干旱区采集大量样本进行样本数据集训练，运用大数据、神经网络、支持向量机进行拟合，通过仿真模拟，选取符合区域特点的参数特征值，提高提取精度，在这方面未来会进行更加深入研究。

参考文献：

[1]周霞，刘彦文，姜宇榕，等. Landsat-8与GF-1遥感影像土地利用数据提取比较——以咸宁市为例[J]. 安徽农业科学，2017，45（31）：213-215，237.

[2]庞新华，朱文泉，潘耀忠，等. 基于高分辨率遥感影像的耕地地块提取方法研究[J]. 测绘科学，2009，34（1）：48-49.

[3]李敏，崔世勇，李成名，等. 面向对象的高分辨率遥感影像信息提取——以耕地提取为例[J]. 遥感信息，2008（6）：63-66，89.

[4]董彦芳，庞勇，许丽娜，等. 高光谱遥感影像与机载 LiDAR 数据融合的地物提取方法研究[J]. 遥感信息，2014，29（6）：73-76，83.

[5]董婷婷，左丽君，张增祥. MODIS数据的水浇地提取[J]. 遥感学报，2009，13（3）：522-534.

[6]张银辉. 农用土地区卫星遥感耕地信息提取及其动态监测技术研究[D]. 泰安：山东农业大学，2002.

[7]白金婷，史园莉. 基于高分影像的面向对象土地利用变化检测方法研究[J]. 农业机械学报，2015，46（增刊1）：297-303.

[8]覃能. 基于高分辨率遥感影像耕地地块提取方法研究[J]. 测绘标准化，2019，35（2）：29-32.

[9]雍万铃. 基于面向对象多尺度分割的目标信息提取研究[D]. 兰州：兰州交通大学，2016：9-16.

[10]徐军. 多尺度分割与案例推理的高分辨率遥感影像信息提取方法[D]. 武汉：武汉大学，2017.

[11]Gang L，Xue X. Fusing texture，edge and line features for smoke recognition[J]. The Institution of Engineering and Technology，2019，13（14）：2805-2812.

[12]Hagara M，Kubinec P. About edge detection in digital images[J]. Radioengineering，2018，27（4）：919-929.

[13]曾建华，黄时杰. 典型图像边缘检测算子的比较与分析[J]. 河北师范大学学报（自然科学版），2020，44（4）：295-301.

[14]Zhou C，Wang G，Huang H，et al. Edge detection based on joint iteration ghost imaging[J]. Optics Express，2019，27（19）：295-307.

[15]张岩，赵玉涛，李光旭，等. 图像的边缘检测技术研究进展[J]. 计算机与数字工程，2020，48（5）：1176-1184，1232.

[16]张紫萍，程小凯，李景文. 基于分水岭和多尺度分割方法的耕地信息提取[J]. 地理空间信息，2012，10（3）：51-53.

[17]崔丽群，张月，田鑫. 融合雙阈值和改进形态学的边缘检测[J]. 计算机工程与应用，2017，53（9）：190-194，200.

[18]王志波，高志海，王琫瑜，等. 基于面向对象方法的沙化土地遥感信息提取技术研究[J]. 遥感技术与应用，2012，27（5）：770-777.

[19]宋浒，张利，许梦晗，等. 基于形态学和小波变换的图像边缘检测方法[J]. 西南大学学报（自然科学版），2020，42（4）：105-111.

[20]吴艳双，张晓丽. 结合多尺度纹理特征的高光谱影像面向对象树种分类[J]. 北京林业大学学报，2020，42（6）：91-101.

[22]严宇宸. 基于纹理特征提取的图像识别方法研究[J]. 电子元器件与信息技术，2019，3（11）：78-79.

[23]高凯旋，焦海明，王新闯. 融合影像纹理、光谱与地形特征的森林冠顶高反演模型[J]. 国土资源遥感，2020，32（3）：63-70.

[24]Canny J. A computational approach to edge detection[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8（6）：79-98.