0 引言

设施农业是现代农业的重要组成部分，作为一种显著延长作物生长季节和提高作物产量的新型生产方式，能够有效提高资源利用率并增强农业系统的防灾减灾能力，是保障“菜篮子工程”顺利开展的重要技术措施^[1]。近年来气候变化导致极端天气事件频发^[2-3]，涝灾、风灾等气象灾害风险的加剧直接威胁设施农业生产^[4-5]。由于建筑结构、材料的差异，不同类型的设施农业对灾害的防御能力完全不同。因此，及时准确地获取设施农业分布情况，可为防灾减灾提供决策支持，并且在促进农业产业结构调整和提高城乡居民的生活水平等方面具有十分重要的现实意义。

随着农业经济、信息技术的迅速发展，设施农业的占地面积逐年攀升。依靠人工实地调查设施农业不仅费时费力，而且很难及时更新大范围设施农业的空间分布数据。与传统地面调查方法相比，遥感技术提供了更高效的手段，并具有成本低、效率高、大范围同步观测的优势^[6-7]，在设施农业的信息提取方面提供了新的思路。目前已有研究尝试使用遥感数据获取设施农业信息，并取得良好进展。提取设施农业通常有2类方法：一是基于像元的支持向量机(support vector machines，SVM)^[8]、随机森林(random forest，RF)和神经网络(neural network，NN)^[9]等机器学习分类方法，其优点是计算效率高且适用性广泛，冯权泷等^[10]利用随机森林模型获取全国塑料大棚分布数据，平均精度为87.45%；二是基于影像分割的面向对象分类方法，例如多尺度分割、灰度分割、知识分割等，可以保证地物在空间上的连续性，但参数依赖性强，仅适合小区域提取，例如吴锦玉等^[11]、Hasituya等^[12]均选择县级区域作为研究区。

现有的研究主要集中在设施农业的整体提取上，但在大范围内针对不同设施农业类型信息提取的方法研究较少。随着遥感技术不断发展，中高分辨率卫星逐步发射，可为地物分类提供更多细节信息。欧洲航天局(European Space Agency，ESA)发布的Sentinel-2光学影像数据空间分辨率达10 m，满足设施农业类型提取需求^[13]。然而影像的空间分辨率越高，同一地物内的异质性越强，从而影响分类精度^[14]。因此除Sentinel-2的光谱特征外，还需考虑其他特征参与精细化分类研究。Sentinel-1搭载的合成孔径雷达(synthetic aperture radar，SAR)可以全天候获取地表信息，提供的后向散射系数能够帮助识别设施农业的纹理特征。融合光学影像和SAR影像，能够实现多源遥感数据互补^[15-16]，进一步提升识别和分类精度。但长时间序列、大空间尺度的遥感数据存储和计算等处理对软硬件环境提出了更高的要求。谷歌于2015年发布了专门用于处理遥感影像的Google Earth Engine(GEE)云计算平台，能够高效管理和处理海量遥感影像，为提取设施农业信息提供了强有力的技术支撑。相比于ERDAS、ENVI等传统遥感图像处理软件，云端处理影像能减少过程数据的生成，用户可以快速实现遥感影像批量处理和结果可视化等功能^[17]，极大地提高了数据处理效率。基于此，国内外学者开展的一系列研究^[18-20]都表明GEE云计算平台在遥感数据分析方面有着巨大的应用潜力。

以山东潍坊为研究区，利用GEE云计算平台，将Sentinel-1雷达影像和Sentinel-2光学影像作为数据源，构建光谱、指数、雷达、纹理和地形特征集合，设计6种设施农业(日光温室、塑料大棚)提取方案，采用随机森林算法分类并进行精度评价，筛选出最佳特征组合，实现大范围提取设施农业类型信息，为设施农业精细化分类提供新思路。此研究可对潍坊设施农业的产业发展和防灾减灾提供重要的决策支持。

1 数据与方法 1.1 研究区概况

潍坊位于山东半岛的西部(35.70°~37.29°N，118.17°~120.01°E；图 1)，属暖温带半湿润大陆性季风气候，降水、光照资源丰富^[21]，适合多种农作物生长。地形以平原为主，兼有丘陵和山地，地势自南向北逐渐降低。南依泰沂山脉，中部为潍北冲积平原，北临渤海莱州湾，地势平坦、土地肥沃。潍坊是中国北方地区最大的蔬菜生产和集散地，设施农业面积位居全省第一。在平原地区聚集不同类型设施农业，分布最为广泛的是日光温室和塑料大棚，设施主要用于种植高附加值的经济作物，如番茄、黄瓜、辣椒、西瓜、草莓等。

1.2 数据 1.2.1 遥感数据获取与预处理

研究使用高空间分辨率的Sentinel-1和Sentinel-2卫星影像，数据源自哥白尼数据空间生态系统(Copernicus Data Space Ecosystem，https://dataspace.copernicus.eu)。Sentinel-1选用干涉宽幅(interferometric wide-swath，IW)模式的地面分辨距离(ground resolving distance，GRD)数据，极化方式选择适用于地表植被特征监测的VV+VH^[22]。Sentinel-2影像数据涵盖从可见光、近红外到短波红外共13个光谱波段，各波段具有不同的空间分辨率(10 m、20 m和60 m)，为提高数据一致性，对所有波段重采样至10 m分辨率。实验收集成像时间范围为2023年1—12月的卫星影像(表 1)。为减少云层遮挡对提取结果的影响，优先筛选云量小于10%的Sentinel-2影像数据，并利用QA60波段进行云掩膜。其中，8月和10月云量过多，无法直接使用单景或多景拼接影像，因此采用中值合成法得到月尺度无云合成影像。

1.2.2 地形数据

研究区地形包括平原、丘陵和山地，加入地形数据的目的是拟提高分类精度。实验基于GEE平台调用美国地质调查局(United States Geological Survey，USGS)提供的30 m空间分辨率的数字高程模型(digital elevation model，DEM)数据产品(USGS/SRTMGL1_003)，重采样至10 m后按研究区边界裁剪，使其与Sentinel数据分辨率匹配，便于多源数据集成与对比分析。

1.2.3 样本采集

在设施农业信息提取过程中，样本类别的划分和选取至关重要。结合实地调查和山东省地理信息公共服务平台天地图2 m历史影像(下载地址为 https://shandong.tianditu.gov.cn)进行目视解译，建立研究区内日光温室、塑料大棚、耕地、道路、水体、森林、建设用地、裸地等主要地物类型的遥感解译标志(表 2)。基于潍坊矢量边界，根据解译标志分散随机绘制2 800个样本，包括日光温室和塑料大棚各580个，其他地物1 640个。按7 ∶3比例划分为训练样本和验证样本，即1 960个训练样本和840个验证样本(图 2)。

1.3 研究方法 1.3.1 技术路线

研究技术路线如图 3所示。使用GEE云计算平台，快速获取覆盖潍坊的Sentinel遥感影像和DEM数据，并进行重采样、去云、裁剪等处理。在此基础上，提取光谱特征、指数特征、地形特征、雷达特征和纹理特征构建特征集。基于不同特征组合，设计6种设施农业提取方案：1.使用光谱特征和指数特征作为基础方案，用以评估其独立的分类性能；2.加入地形特征，提供额外的空间位置信息，有助于提高分类模型的判别能力；3.在基础方案中加入雷达特征用于反映地物表面粗糙度、电磁波散射特性等属性，是对光学影像的重要补充；4.引入基于Sentinel-2影像计算的纹理特征，捕捉地物结构的细微差异，从而提升设施农业精细化分类能力；5.综合全部特征，全面利用多源数据的优势最大化提升分类精度；6.在全部特征组合的基础上，通过特征优选方法筛选贡献度高的变量，剔除冗余信息，进一步提升分类效率和精度。所有方案均采用随机森林监督分类方法进行提取，并对分类结果进行精度评价，选择精度最高的特征组合方案，获取潍坊设施农业空间分布信息。

1.3.2 特征选取

(1) 光谱特征。在Sentinel-2光学影像的13个原始波段中，剔除B1、B10这2个无关波段，选取其他11个光谱波段作为原始波段特征，并计算19个常用遥感指数(表 3)，如归一化植被指数(NDVI)、归一化水体指数(NDWI)、红边位置指数(REPI)、土壤调节植被指数(SAVI)等。计算结果在实验中作为独立的指数特征波段，与原始波段并列参与分类任务。

(2) 雷达特征。Sentinel-1有4种成像模式(干涉宽幅模式IW、扩展宽幅模式EW、波束模式WV、条带模式SM)和2类极化方式，分为单极化方式(HH或VV)和双极化方式(HH+HV或VV+VH)。其中，IW模式广泛应用于陆地表面的遥感研究^[23]，特别是VV和VH偏振组合，在地表特征提取和分析中表现出优异的效果。本研究选用VV+VH作为独立的雷达特征波段参与特征构建。

(3) 纹理特征。纹理特征是地物提取和识别中的重要参数，用来描述图像中像素灰度值的空间分布，这种局部的统计信息能够帮助提取复杂地物的特征。常用的纹理分析方法是灰度共生矩阵(GLCM)^[24]，因其计算过程简单、结构信息丰富，故在遥感提取研究中得到广泛应用。GLCM算法需要8位灰度图像作为输入，因此此实验通过原始影像的红(RED)、绿(GREEN)、蓝(BLUE)波段的线性组合生成所需灰度图像(gray)，计算公式如下。

$ \text { gray }=(0.3 \times \text { RED })+(0.59 \times \text { GREEN })+(0.11 \times \text { BLUE }) $

(1)

实验选取8种主要的纹理特征，包括角二阶矩(ASM)、相关性(Correlation)、对比度(Contrast)、熵(Entropy)、均值(Mean)、方差(Variance)、相异性(Dissimilarity)、同质性(Homogeneity)。

(4) 地形特征。日光温室和塑料大棚的建造受地形限制，加入地形特征尝试提高分类精度。除海拔高度(Elevation)外，增加3种常用的地形特征，即坡度(Slope)、坡向(Aspect)、山体阴影(Hillshade)，这4种特征将作为独立波段参与特征构建。

1.3.3 Jeffries-Matusita距离

Jeffries-Matusita(JM)距离^[25]常用于遥感影像分类中评估不同类别之间的可分性。将全年不同月份的JM距离进行比较，JM距离越大，表明设施农业与其他地物类别的可分性越好，在此研究中JM距离最大的月份被认为是提取的最佳时相。

设μ_i和μ_j分别为类别i和j的均值向量，Σi和Σj分别为对应的协方差矩阵。JM距离的定义如下。

$ \mathrm{JM}_{i j}=\sqrt{2\left\{1-\exp \left[-\frac{1}{8}\left(\boldsymbol{\mu}_i-\boldsymbol{\mu}_j\right)^{\mathrm{T}}\left(\frac{\sum i+\sum j}{2}\right)^{-1}\left(\boldsymbol{\mu}_i-\boldsymbol{\mu}_j\right)-\frac{1}{2} \ln \left(\frac{\frac{\sum i+\sum j}{2}}{\sqrt{\left|\sum i\right|\left|\sum j\right|}}\right)\right]\right\}} $

(2)

其中：μ_i－ μ_j表示类别均值向量的差异，$\frac{\sum i+\sum j}{2}$表示类别协方差矩阵的均值。

1.3.4 累计差

在计算纹理特征时，滑动窗口大小(邻域的尺寸)对分类结果有显著影响。较小窗口容易引入过多噪声，导致分类结果不稳定。而较大窗口会模糊图像中的细节特征，降低分类精度。因此，选择合适的窗口大小是获取纹理特征的第一步。此研究采用累计差(Δf)^[26]衡量地物区分度以确定窗口大小，累计差越大，表明地物在当前窗口下的区分度越高，其计算公式如下。

$ \Delta f=\sum\limits_{i=1}^k\left(\left|u_i-u_{i+1}\right|+\left|u_i-u_{i+2}\right|+\cdots \cdots+\left|u_i-u_{i+k-1}\right|\right) $

(3)

其中：Δf代表某一纹理特征变量的累计差，此研究共有8个纹理特征，因此f=1，2，……，8；k为地物类别总数；u_i表示第i类地物在某特征变量上的均值。

1.3.5 分类方法

RF算法^[27]是一种基于集成学习的监督分类算法，其核心思想是通过构建大量相互独立的决策树来提高分类的准确性和稳定性，具有良好防止过拟合能力。相比其他分类器，RF在处理高维遥感影像数据时表现出更好的性能，对噪声和缺失数据具有较强的鲁棒性，适用于多类别和多特征的遥感分类^[28]。

在GEE云计算平台中，使用RF分类算法时需要设置关键参数，此研究设置决策树的数量为100个，每个决策节点使用的变量数为特征总数的平方根。

1.3.6 精度评价

在遥感图像分类中，混淆矩阵是评估分类模型性能的常用工具，用于表示实际类别与模型分类结果的对比情况。此研究通过计算分类结果的混淆矩阵对6种特征组合方案进行精度评价，评价指标包括总体精度(overall accuracy，OA)、Kappa系数、用户精度(user’s accuracy，UA)和生产者精度(producer’s accuracy，PA)。

2 结果与分析 2.1 最佳时像选取

根据目视解译的样本，使用GEE云计算平台，对Sentinel-2影像不同月份下设施农业与其他地物的JM距离进行计算，确定提取的最佳时相。如图 4所示，2023年Sentinel-2影像塑料大棚与其他地物的JM距离均大于日光温室与其他地物的JM距离，且二者在全年变化趋势相似，夏季JM距离小，冬季JM距离大，2月20日均达到最大。这表明设施农业在2月的光谱差异性最大，更容易与其他地物区分。因此，选取2月20日的Sentinel影像进行设施农业提取。

2.2 纹理特征提取

GEE云计算平台提供的GLCM算法中可定义计算窗口大小。研究选取3×3、7×7、11×11、15×15、19×19、23×23、27×27、31×31、35×35、39×39、43×43、47×47、51×51、55×55共14种尺寸的滑动窗口，在最佳提取时相的Sentinel-2灰度图像上计算不同窗口大小下各地物(日光温室、塑料大棚及其他地物)的纹理特征累计差。如图 5a所示，累计差随着窗口的增大，先增加后减小，在27×27窗口下累计差达最大，此时设施农业与其他地物的纹理特征区分度最高，选择27×27窗口为最佳提取窗口。图 5b展示了最佳提取窗口下各纹理特征的累计差，其中方差特征的累计差最高，熵、相异性、对比度、均值依次降低，同质性、角二阶矩、相关性的累计差最低且与前5个纹理变量相差过大，累计差均小于0.5。因此选择方差、熵、相异性、对比度、均值作为纹理特征参与模型训练。

2.3 特征对分类的贡献分析

为了筛选出对分类精度贡献较大的特征，避免不重要的噪声特征干扰，使用基于基尼不纯度的特征重要性评估，实验共获取特征变量41个。由图 6可以直观看出，地形特征中的Elevation重要性最高，这反映出设施农业的分布可能与地形起伏密切相关，日光温室和塑料大棚多修建于平坦的地理条件下，与潍坊的实际情况相符合。后面依次是均值、方差、相异性、熵，表明纹理特征对此次分类重要性很高，这对于捕捉设施农业的空间结构具有重要作用。此外，雷达波段VV、VH和光谱波段B9、B12也具有较高的重要性。光谱特征尤其是与近红外和可见光相关的波段常用于植被监测，在提取与设施农业相关的植被信息时起到了关键作用。特征重要性排序依次为：纹理特征、光谱特征、雷达特征、指数特征、地形特征。

特征冗余不仅会增加计算复杂度，还会出现过拟合的风险，导致分类精度降低，因此需要对特征进行筛选。根据重要性顺序依次增加特征训练RF模型，计算每个模型分类结果的OA。由图 6可以看到，添加NDVI后，整体分类精度趋于平稳且无显著提升，因此保留包括NDVI在内的前26个特征构成最优特征集。

2.4 精度评价

使用固定种子的RF算法对6种特征组合方案进行分类，结果如表 4所示。方案1仅包括光谱和指数特征，分类精度在各方案中最低。尽管光谱信息和植被指数是常用的遥感特征，但对于复杂的设施农业，并不足以完全反映其特征差异。在方案1基础上分别加入地形特征、雷达特征、纹理特征后，分类精度均有提高。其中，雷达特征对分类精度的提升最少，塑料大棚的UA提高2.28%，而日光温室仅提高0.54%，说明雷达特征对日光温室提取精度的贡献有限。在方案1中增加纹理特征后，OA提升至88.78%，且日光温室和塑料大棚的分类精度均有大幅提升。对比方案2和方案3，纹理特征对分类精度的提高最显著，表明纹理特征能够捕捉影像中局部空间关系和细节，对于区分形态复杂的设施农业尤其有效。方案5融合所有特征，OA达89.38%，经过优选特征后，方案6的OA相较方案1提高了6.33%，Kappa系数达0.86，有效减少特征冗余，进一步提高了模型的鲁棒性和分类精度，在所有方案中取得最好的分类效果。根据OA和Kappa系数，不同特征对设施农业提取的有利程度排序为：纹理特征、地形特征、雷达特征、光谱特征。

基于GEE云计算平台和随机森林算法，提取出潍坊2023年设施农业的分布区域。图 7展示了日光温室和塑料大棚集中分布区域、交错分布区域、复杂背景区域等4个典型区域的提取情况。对于单一类型且大面积分布的设施农业，所有方案的分类结果差异较小。然而，在交错分布和复杂背景区域，6种方案表现出较大差异。方案1—3均可识别设施农业的基本轮廓，但在边缘处出现明显的分类错误和噪点，导致边界不清晰、误分类现象较多。方案4纹理特征的加入显著提高了分类精度，特别是在细节处，混合像元更多的被分为其他地物，使道路、农田的边界更加清晰，但会使设施农业统计面积减少。方案5综合了所有特征，包括光谱特征、指数特征、地形特征、雷达特征和纹理特征。虽然该方案理论上应达到最优效果，但实际分类结果并未明显优于方案4。这可能是特征过多导致信息冗余，引入了额外的噪声。而特征优选后的方案6在不同区域下均能够较好地捕捉设施农业的完整性和边界特征。

图 8展示了特征优选方案的设施农业提取结果。可以看出，潍坊的设施农业主要集中在寿光和青州，总面积占全市设施农业面积的66%(表 5)，但日光温室和塑料大棚的空间分布存在明显差异。日光温室主要集中分布在寿光、青州、临朐等区域，其中寿光面积达40.36万亩^①，占全市日光温室总面积的51.45%，是日光温室设施农业发展的核心区域。其次为青州，面积为17.09万亩，占21.79%。与日光温室相比，塑料大棚的分布更分散。青州和昌乐的塑料大棚面积较大，分别为16.89万亩和8.79万亩，合计占全市塑料大棚总面积的57.52%。

① 1亩≈0.066 7 hm²。

与地形图叠加发现，设施农业主要建在地势平坦的平原地区，极少建立在山地和丘陵地区。从地形地貌的角度来看，寿光南部和青州东北部地势相对平坦，主要为冲积平原，土地开阔平整，减少了建设温室大棚基础设施的成本，有利于大规模设施农业的布局和管理^[29]。在人文因素方面，寿光作为“中国蔬菜之乡”，拥有丰富的设施农业经验和完善的产业链条。青州北部和昌乐北部由于与寿光接壤，很大程度上受其技术和市场的辐射和影响，形成了设施农业集聚区。

3 讨论

借助GEE云计算平台，基于Sentinel-1/2和DEM数据，利用RF算法，设计6种特征组合方案，实现对潍坊设施农业信息的提取。与单一数据源研究^[10]相比，此研究充分利用了多源遥感数据的融合优势。朱德海等^[30]、汤紫霞等^[31]仅基于光学卫星影像结合纹理特征对大棚进行提取，而此研究在此基础上加入雷达数据，使平均分类精度显著提高。目前关于设施农业分类提取的研究较少，且均采用面向对象的分类方法。例如王乌云等^[32]利用面向对象多尺度分割方法对内蒙古河套灌区大棚信息进行提取。但面向对象的分类高度依赖于影像的分割质量，对空间分辨率、研究区范围较为敏感。识别精度受设施类型、分布情况和面积占比不同的影响^[33]。此研究采用RF算法结合特征优选，在减少计算复杂度的同时，提高泛化能力，保证高分类精度。这些改进为设施农业的遥感提取提供了高效、可靠的技术方案，并为设施农业灾害风险管理提供决策支撑。

日光温室与塑料大棚的分类精度易受“同物异谱”和“异物同谱”现象干扰。两者通常采用聚乙烯薄膜作为覆盖材料，在Sentinel-2影像的可见光-近红外波段反射率差异较小，尤其在薄膜老化或污损后可分性进一步降低。光谱特征的局限性促使研究引入其他特征以增强分类能力。加入纹理特征后，日光温室和塑料大棚的PA和UA与基础方案1相比分别提升7.32%和7.02%、5.29%和5.81%，与塑料大棚相比，日光温室提取效果受纹理特征影响更大，主要考虑2方面原因：(1)日光温室的建筑结构较为复杂，包含厚实的背墙和两侧墙体、透明材料的屋顶斜面、多层的保温棉被等^[34-35]，这些具有不同光照反射和吸收的地物都会在遥感影像中形成明显的纹理特征。而塑料大棚结构相对简单，表面有较为均匀的塑料薄膜覆盖^[35]，纹理变化较少。因此，纹理特征对日光温室的捕捉能力更强，能够更清晰地区分日光温室与其他地物，使得分类精度提高。(2)在建造日光温室时，通常在温室之间留出室外作业通道和灌排水渠道^[36]，在遥感影像上呈现出较复杂的背景，纹理特征更加明显。而塑料大棚通常成片布局，表面有相对统一的光滑程度，受周围背景的干扰较小，因此纹理特征对塑料大棚的识别贡献较小。

此外，混合像元效应也是影响分类精度的重要因素。当地物边界与像元网格未对齐时，其边缘像元将包含设施农业与背景地物(如裸土、道路)的混合光谱信号，导致模型无法有效分离地物特征(图 7)。这一现象可能导致尺寸较小的单体大棚不能被准确提取或漏提。未来研究可以考虑融合亚米级卫星数据提高空间分辨率，或构建多尺度特征学习模型以抑制混合光谱干扰。

尽管特征优选方法在设施农业类型的提取中取得了较高的精度，但研究中仍有不足的地方。首先，此研究主要集中在单一时相影像的提取，未充分考虑多时相数据对提取精度的影响，进一步的研究可以结合设施农业内部种植作物的种类，加入多时相数据，更好地捕捉季节性变化和设施农业使用状态。其次，随着人工智能的不断发展，深度学习技术在许多领域表现出优异的性能，未来研究可以探索利用深度学习算法，实现大面积、更加精细的设施农业信息提取。

4 小结

(1) 根据样本计算不同设施农业类型与其他地物的JM距离，得出日光温室和塑料大棚的最佳提取时相均为2月。

(2) 使用累计差法从14种窗口大小中确定提取设施农业的最佳纹理特征窗口是27×27，最佳纹理特征组合为均值、熵、方差、相异性和对比度。

(3) 不同的特征组合方案对提取结果影响不同。分类特征对设施农业的重要性排序依次为：纹理特征、地形特征、雷达特征。3个特征的加入使得分类精度分别提升4.81%、1.66%、0.94%。

(4) 6个方案中特征优选方案提取不同设施农业类型信息效果最佳，OA为91.03%，Kappa系数为0.86。其中，塑料大棚的PA和UA分别为92.21%和91.84%；日光温室的PA和UA分别为88.54%和89.47%。这表明该方法能够有效获取不同设施农业类型的空间分布和面积信息。