遥感已然成为土地监测的重要手段^[1-4]，利用遥感影像进行人工种植园提取对于农作物估产、林业调查、土地评估等具有十分重要的意义^[5]。人工种植园的提取主要依据有两点：一是物候学特征，不同植被类型在不同的生长发育阶段表现出迥然不同的形态；二是利用不同植被类型表现出的细微光谱反射差异。因此现阶段人工种植园的提取主要依靠多时相数据和高光谱数据。由于不同类型的植被往往难以用宽波段遥感影像分开，因此对于宽波段遥感影像数据，种植园的提取需要融合光谱外更多的信息。文献[6]基于面向对象思想，对各分割对象计算52个特征，其中包括光谱、纹理以及多种地形辅助数据，并采用历史航空影像解译结果、GPS辅助数据以及野外调查目视解译结果3种方法获取样本进行分层分类，初步实现了植被的精细化分类。由于自然植被对于环境的依赖性，因此地形环境的辅助数据在分类中起到了重要作用。目前，种植园的提取对数据源具有较强的依赖，在特征上主要采用以植被指数为代表的光谱特征以及灰度共生矩阵等微纹理，对于其排列结构的规律性模式并未得到足够的重视。

结构纹理作为大尺度宏纹理，通常用来提取纹理基元及其排列模式。例如，文献[7]基于灰度差分统计来检测各向异性、对称性、规律性。文献[8]基于二值图像计算多方向和尺度的灰度共生矩阵，利用其中几个特征的最大值来提取网格结构的位移矢量。文献[9]基于自相关函数提取规则纹理的纹理基元和周期。文献[10]利用傅里叶谱的两个基频来鉴别形成网格结构的纹理基元。这些结构方法要求纹理单一且具有严格的排列规律，对于植被等复杂多变的自然地物适用性较弱。文献[11]通过对单株树木的检测，计算树冠半径、邻域树距离等格网测量模式及树冠色调饱和度、NDVI均值方差等多维特征，融合多分类器实现了柑橘园85%自动分类和15%目视解译。但是该方法对于有杂草和高度植被化的土地背景，单株树木的检测及其与背景的分离十分困难。

本文研究目的是探索种植园的纹理结构，提出一种结构指数特征，能够满足多尺度多方向种植模式，允许纹理基元和格网种植线的变形，降低对数据源的要求，并且具有高分离能力，从而减少特征冗余，简化提取过程，实现依靠单一或者少数特征的种植园自动提取。

1 种植园特性分析

光谱特征、时相特征和空间特征是遥感地物提取的3大基本特征。与其他绿色植被一样，种植园具有典型的植被光谱反射特性曲线，在蓝光和红光波段有两个吸收带，在绿光波段有一个明显的反射峰，至近红外波段范围内反射率达到高峰。因此，近红外影像常常用来进行植被提取。而对于种植园与其他植被的区分，单纯依靠宽波段遥感影像的光谱特征往往难以取得理想效果，需借助高光谱影像对不同植被类型的光谱细微差异进行更深入的分析。种植园的时相特征主要表现为种植园一般具有特定的生长规律和物候特征，因此可以利用多时相影像数据进行种植园提取。但是，无论高光谱影像还是多时相影像，数据获取成本均较高。其中，高光谱影像往往空间分辨率不高且易出现“维数灾难”现象^[12]，多时相影像重访周期长，受云雨天气影响，往往难以获得“最佳物候期”的图像数据^[13]，且多时相遥感影像必须进行精确的几何配准，研究表明，即便是微小的配准误差，也会大大降低结果的精度^[14-15]。此外，对于影像中大量存在的同物异谱和异物同谱现象，单纯依靠光谱特征和时相特征往往难以有效解决。种植园空间特征的有效提取能够极大地抑制同物异谱和异物同谱问题，降低对数据的依赖，以同时满足低成本、高时效的种植园提取。

株距和行距是普遍的种植园空间布局模式。无论什么类型的种植园，为了保证种植园CO₂浓度、温度、湿度、太阳有效辐射和净辐射等以利于种植物生长，种植园的布局往往是保留一定的株距和行距，且同一片种植园的株距和行距一般较固定。与其他地物相比，这种近乎等间隔的空间布局构成种植园独有的标志性特征。在此理论基础上，本文提出一种多尺度多方向结构指数对种植园空间布局进行数学描述，通过构造多尺度多方向窗口垂直投影，获取地物窗口投影曲线，种植园由于具有近乎等间隔的排列结构，因而其投影曲线具有较规律的波形，波峰和波谷分别对应植株和株距中心，二者在横轴上的位置间隔较稳定，因此通过相邻波峰波谷之间距离的归一化二次差分来进行度量，并结合形态学后处理实现种植园提取。

2 多尺度多方向结构指数及种植园自动提取算法 2.1 LOG滤波预处理

为了凸显种植园纹理基元，需要对影像进行预处理。由于单株种植物即纹理基元一般呈中心向四周发散状，例如柑橘树的树冠在遥感影像上类似“圆形”，从中心散开。因此选择对各个方向的反差均能敏感检测的各向同性滤波器，高斯拉普拉斯函数(Laplace of Gaussian function, LOG)^[16]，进行预处理。LOG算子源于文献[16]中提出的边缘提取思想，即首先对图像进行最佳平滑，最大限度地抑制噪声，然后求取边缘。其定义如下

式中，Δ²为拉普拉斯算子；G(x, y)为高斯函数；σ为高斯分布因子。由式(3)可见，σ的取值变化对x和y的影响相同，体现出圆对称的性质。σ值越大，滤波器对图像的平滑作用越明显。本文选择σ值分别为5、7、9，对图像进行多尺度滤波，以满足种植基元不同尺度的需求。

2.2 窗口纹理结构分析

假设种植园方向与窗口方向一致，如图 1(a)所示，LOG滤波后的窗口影像如图 1(b)所示。窗口纹理结构分析步骤如下：

(1) 通过式(4)将窗口内所有像素按列沿窗口垂直方向进行累加，获取窗口投影曲线^[17]，如图 1(c)所示，种植园波形较规律，而非种植园区域波形紊乱。式(4)计算公式如下

式中，x_c为第c列的垂直投影值；width为窗口宽度；G_r为第c列第r行的灰度值；height为窗口高度。

(2) 窗口垂直投影曲线上，种植园区域的波峰和波谷分别对应单株种植物中心和种植间距中心，为了准确获取二者在横轴上对应位置和其距离值，通过计算局部极值对窗口垂直投影曲线进行波谷波峰分割，如图 1(d)所示，并获取相邻波谷波峰之间的距离dist_i(i=1, 2, …, n)。

(3) 通过式(5)计算曲线段规律性值y_i。式(5)中，为了对等间距的种植园结构进行准确度量，y′_i取相邻4对波谷波峰之间距离dist_i的二次差分，并进行归一化，如式(6)和式(7)所示，对于完全等间距的种植园区域，y′_i的值为1，而由于各种因素，实际间距往往近似相等，经过归一化处理，相邻距离二次差分层层缩小，最后获取的y′_i值也接近1。非种植园区的曲线往往会出现相邻波峰波谷之间的距离和投影幅值差过大及过小的情况，前者与纹理基元尺度和种植间隔相关，考虑种植紧密的小尺度区域，下限不作约束，而后者与窗口大小相关，上限不作约束，因此，分别设置相邻波峰波谷距离差和幅值差的约束参数ρ、ε，滤除部分非种植园区域。最后通过上采样得到每一个曲线段对应的规律性值y_i∈[0,1](i=1, 2, …, n)，如图 1(e)所示，种植园区域和非种植园区域y_i值分别逼近1和0。式(5)-式(7)计算公式如下

(4) 将曲线段规律性值y_i(i=1, 2, …, n)映射到对应的窗口影像像素特征r，获取结构特征图，如图 1(f)所示，除了少量孤立片段，对应的种植园区域与非种植园区域反差明显。通过形态学开运算进一步剔除小面积孤立片段，最终获取窗口影像结构特征，如图 1(g)所示。

2.3 多尺度多方向结构指数

为检测种植园方向，并且考虑多尺度多方向种植园区，采用多尺度多方向投影窗口进行纹理结构分析，具体步骤如下：

(1) 对影像进行多尺度LOG滤波，如2.1节所述，LOG滤波器中，高斯滤波尺度分别选取5×5、7×7、9×9，拉普拉斯滤波尺度固定为3×3。

(2) 构造多方向搜索窗口，如图 2所示，窗口方向θ取0~180°，步长间隔1°，窗口高度height根据影像尺度自定义，过大容易造成误提取，过小容易造成漏提取。

(3) 在各尺度各方向下，按照2.2节所述方法，逐一进行窗口纹理结构分析，获取对应的像素特征值r_{[scale][angle]}, (scale∈[5,7,9], angle∈[0,180])。

(4) 利用式(8)对各像素特征值r_{[scale][angle]}进行最大值滤波，获取多尺度多方向结构指数(multi-scale and multi-direction structure index，MMI)

2.4 种植园自动提取

在MMI特征图上，种植园区域和非种植园区域的特征值分别逼近1和0，二者反差显著。良好的可区分性使得MMI能够凭借简单的阈值实现对种植园的提取。具体算法为，通过对MMI设固定阈值0.8，即当像素的MMI特征值大于0.8时，对种植园进行标记提取。然后通过形态学后处理对孤立小面积的误标记区域进行剔除。形态学腐蚀能够快速有效剔除小面积非种植园区，但是同时会破坏提取的种植园区，因此在进行形态学腐蚀之后，利用原提取图进行形态学重建，如式(9)所示，实现在保留正确提取基础上的误提取剔除

式中，R表示形态学重建；e表示形态学腐蚀；s为结构算子的大小；MMI为原提取特征图。

3 试验 3.1 试验设计

试验部分通过对3组不同种植模式、背景植被化程度、结构复杂度的种植园进行试验以分析本文算法的有效性和适用性。试验数据为6幅分辨率0.5～0.65 m的多光谱影像，为了验证本文算法对数据的依赖性以及光谱相关性，将其中3幅多光谱影像去彩色处理获取单通道灰度影像，每组试验分别利用两种影像，对本文算法和传统算法进行分析比较。种植园提取结果用绿色进行标记。地面参考影像为分别对试验影像中种植园(绿色标记)和非种植园(白色标记)进行随机密集采样所得，以统计种植园提取结果的误分和漏分情况，漏分率为试验结果对于采样像素的漏分数量与总体采样数量的比值，误分率为试验结果中误分的像素数量与提取的植被像素总数的比值。

3.2 试验结果与分析

(1)第1组试验中多光谱影像(图 3(a))和灰度影像(图 3(d))，分辨率0.5 m，图中种植园总体结构呈近似网格排列，种植线存在多方向，种植基元存在多尺度且具有明显变形。多光谱影像的试验参数设置为height=26、ρ=125、ε=6，灰度影像为height=26、ρ=125、ε=6。

针对网格结构的种植园提取，传统的方法存在的局限是假设图中只存在单一纹理，且具有理想的纹理基元大小和严格的网格排列结构^[18-20]。而第1组试验中存在多种类型种植物，呈多方向排列，且单株种植物的大小和植株间隔不固定，因此其纹理结构具有多尺度多方向性，纹理基元存在明显变形，网格结构不规律。由图 3(b)和图 3(e)可见，本文方法突破了传统方法的局限，对其能够进行有效提取，分别基于地面参考影像图 3(c)和图 3(f)进行精度计算，多光谱影像提取结果漏分率1.91%，误分率1.85%，整体精度98.12%，灰度影像提取结果漏分率2.34%，误分率0.82%，整体精度98.42%。

(2)第2组试验中多光谱影像(图 4(a))和灰度影像(图 4(d))，分辨率0.5 m，种植园成近似直线排列，种植紧密，背景高度植被化。多光谱影像试验参数设置为height=30、ρ=130、ε=6，灰度影像参数设置为height=30、ρ=150、ε=6。

由于种植紧密且背景高度植被化，单株树木与背景分离困难，加上单株种植物的长势各异，因此传统基于单株种植物检测和测量的种植园提取方法^[11]并不适用。由图 4(b)和4(e)可见，无论是多光谱影像还是仅仅具有单通道的灰度影像，本文方法对此能够有效提取，分别基于图 4(c)和4(f)进行精度评价，多光谱影像提取结果漏分率3.89%，误分率1.04%，整体精度97.54%，灰度影像漏分率3.85%，误分率1.15%，整体精度97.5%。

(3) 第3组试验中多光谱影像(图 5(a))和灰度影像(图 5(d))，分辨率0.6 m。两幅影像中种植园规律性较差，空间结构无固定方向和形态，并且种植园比较破碎，背景高度植被化。多光谱影像的试验参数设置为height=30、ρ=130、ε=6，灰度影像为height=30、ρ=130、ε=6。

这是极为复杂的一类种植园，从光谱上与背景分离困难，单株植被难以被检测，从纹理结构上，规律性差，方向和形态多变，缺乏整体性固定的排列模式，且较为破碎，图 5(a)中种植园更是呈零散分布，对种植园提取算法是一种相当大的挑战。目前，对于这种种植园的提取，多是依靠人工目视解译。本文提出的MMI能够对其复杂的空间结构进行有效表征，并且作为唯一特征实现种植园自动提取，如图 5(b)和图 5(e)所示。根据随机密集采样获取的地面参考影像(图 5(c)和图 5(f)), 种植园提取结果精度评价分别为：多光谱影像提取结果漏分率15.28%，误分率0.35%，整体精度92.19%，灰度影像提取结果漏分率10.74%，误分率2.33%，整体精度93.47%，显示了MMI对种植园纹理结构特征描述的有效性和本文方法具有较强的适用性，对于极复杂的种植园也能进行提取，并达到90%以上整体精度。

4 结论

本文提出了一种能够有效描述种植园纹理结构的特征，多尺度多方向结构指数(MMI)，并且结合形态学后处理实现了种植园的自动提取。本文算法简单、自动化程度高，3组试验灰度影像提取结果精度与多光谱影像相当，说明MMI特征光谱相关性低，且仅依靠MMI单一特征实现提取，无需其他特征和辅助数据，从而降低了对数据的依赖。通过设计各种类型种植园提取试验显示，本文方法突破了传统方法的限制，具有较强的适用性，对于多种复杂的种植园区域均能实现有效提取，总体精度达到90%以上。

本文的不足之处在于试验参数的设置，包括窗口高度height，相邻波峰波谷之间的距离上限ε和投影幅值差下限ρ，以及形态学后处理中结构算子的设置，均为经验性值，与影像分辨率、种植园尺度等有关，但尚无可以计算的依据，需要进一步研究改进。