吉安市地处江西省中西部, 位于赣江中游, 主要以山地和丘陵地貌为主, 如图1(a— e)所示, 境内四季分明、 雨量充沛, 为亚热带季风气候, 年日照时数为1 600~1 800 h, 年降水量在1 400~1 600 mm之间, 年平均温度约17 ℃, 极适宜油茶树生长。 在中央财政与当地政府的鼓励下, 吉安市的油茶产业近年来得到了快速发展, 以油茶为主要经营树种的扶贫造林专业合作社超过10家, 油茶企业信息备案登记在册有19家, 油茶林已成为当地农民的重要收入来源之一。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区示意图Fig.1 Location of the study area

选用Sentinel-1 SAR雷达遥感数据和Sentinel-2A光学遥感数据提取油茶林种植区域。 Sentinel-1 SAR是主动微波遥感卫星, 提供四种极化模式(VV、 VH、 HH、 HV), 不受天气、 光照等条件的影响, 对云层具有较强的穿透力。 Sentinel-2A是多光谱遥感卫星, 其搭载了一台多光谱成像仪, 可在13个光谱波段范围内获取图像数据, 捕捉地表的细微特征。

Sentinel-1 SAR雷达遥感数据对大气层的穿透力强, 对地面纹理特征敏感。 由于其成像方法和特点, 目标研究区域的地形会对成像产生干扰和噪声, 不利于后续油茶林分类。 本研究采用Lee滤波器对Sentinel-1 SAR中所用到的波段进行滤波。 Lee滤波器通过计算局部窗口内的平均值, 并结合每个像素点的权重系数, 实现对噪声的抑制。

Sentinel-2A光学遥感影像在进行农作物时空分布信息提取时, 由于数据量化、 传感器限制等因素, 导致图像中地物特征和边界的细节不够突出, 影响油茶林遥感识别和分类。 本研究先对Sentinel-2A光学遥感影像进行去云处理, 再采用基于梯度卷积的Sobel边缘检测算子对影像进行边缘检测处理以增强影像中的高频细节。 Sobel算子通过计算图像中每个像素点周围像素的灰度加权差检测边缘。 图2(a)为遥感影像, 图2(b)为Sobel边缘检测算子效果, 遥感影像用于对比。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 遥感影像(a)与Sobel边缘检测算子效果图(b)Fig.2 Remote sensing image with Sobel edge detection operator effects

采用谷歌地球专业版云平台提供的点取样方法和实地调查。 对油茶林样地进行了实地调查。 其余样地结合谷歌地球专业版在线选取并验证样本准确性, 然后导入 GEE 云平台。 在线选择过程基于多样性、 随机均匀性和可验证性原则, 以确保研究区域内样本的空间分布尽可能满足代表性和均匀性的特点。 为准确反映模型效能, 从中随机抽取70%作为训练样本, 30%作为测试样本。 研究中五种土地类型样本点信息见表1。

表1

Table 1

表1(Table 1)

表1 样本数据集 Table 1 Sample data sets 数据类型 水体 不透 水面 其他 林地 耕地 油茶 林 合计 训练数据集数量 439 654 970 750 339 3 152 测试数据集数量 188 281 416 322 146 1 353 总计 627 935 1 386 1 072 485

表1 样本数据集 Table 1 Sample data sets

本研究技术路线图如图3所示。 首先采用GEE云平台2021年的Sentinel-1 SAR和Sentinel-2A遥感影像数据, 构建地物类型的光谱特征、 纹理特征、 地形特征和极化特征, 采用Sigmoid激活函数对光谱特征进行处理, 采用多种机器学习分类器和卷积神经网络对2021年油茶林样本点进行训练和精度验证, 并通过基于余弦相似度的样本迁移, 分析2019年至2023年的吉安市油茶林种植面积变化趋势及空间分布。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 技术路线图Fig.3 Flow chart of proposed method

1.4.1 特征构建

使用Sentinel-1 SAR卫星的VV和VH波段构建极化特征能够提供关于地表覆盖和物理特性的重要信息。 其中VV极化表示垂直发射和垂直接收的极化方式, 主要用于表征地物散射的垂直分量, 对于具有垂直极化特征的地物具有较高的敏感度。 VH极化表示垂直发射和水平接收的极化方式, 主要用于探测地物的极化散射特性和地表的粗糙度, 能够提供更多关于地表物理特性的信息。

结合研究区的地理位置与油茶的物候期, 选取 Sentinel-2A卫星搭载的多光谱传感器获得多光谱波段数据。 通过对光谱波段数据进行数学运算或组合得到相应的指数特征, 进而构建特征集如表2所示。 归一化差值建筑指数(normalized difference built-up index, NDBBI)通过建筑材料和裸土在短波红外和近红外波段的反射差异来识别建筑区域, 用于城市扩展监测和土地利用变化检测。 土壤水分指数(soil water index, SWI)通过短波红外和近红外波段的反射差异来估算土壤水分含量, 用于农业和水资源管理中的土壤湿度监测。 大气阻抗植被指数(atmospherically resistant vegetation index, ARVI)通过使用蓝光波段对红光波段进行校正, 减少大气气溶胶的影响, 用于植被监测, 特别是在大气条件不稳定的地区。 增强型建筑与裸地指数(enhanced built-up and bareness index, EBBI)通过短波红外和红光波段的反射区分建筑物和裸土, 用于城市扩展监测和裸土识别。 水体提取指数(water extraction index, WEI)通过绿光和近红外波段的反射差异来识别水体, 用于水体监测和水资源管理。 增强型植被指数(enhanced vegetation index, EVI)增强了植被信息, 特别是在高生物量区域。 相比NDVI, EVI对大气条件和土壤背景具有更好的校正能力, 用于植被类型识别和覆盖度估算。

表2

Table 2

表2(Table 2)

表2 光谱指数特征公式 Table 2 Spectral index characteristic formula 指数 公式 所用波段 NDBBI NDBBI= 1.5×ρSWIR2-(ρNIR+ρGREEN)÷21.5×ρSWIR2+(ρNIR+ρGREEN)÷2 B3, B8, B12 SWI SWI= (ρREDEDGE1-ρSWIR1)(ρREDEDGE1+ρSWIR1) B5, B11 ARVI ARVI= ρNIR-(ρRED-0.5×(ρRED-ρBLUE))ρNIR+(ρRED-0.5×(ρRED-ρBLUE)) B2, B4, B8 EBBI EBBI=0.08× ρ NIR+0.52× ρ SWIR1+1-0.25× ρ GREEN-0.02× ρ RED B3, B4, B8, B11 WEI WEI= (ρGREEN-ρNIR-ρSWIR1)(ρGREEN+ρNIR+ρSWIR1) B3, B8, B11 EVI EVI= (ρNIR-ρRED)(ρNIR+6×ρRED-7.5×ρBLUE+1) B2, B4, B8

表2 光谱指数特征公式 Table 2 Spectral index characteristic formula

有研究表明农作物遥感分类研究仅计算光谱特征会导致“ 同物异谱、 异物同谱” 的现象, 影响最终分类效果, 而不同地物类型在遥感影像上纹理表现差异较大。 油茶林在遥感影像上具有规则的纹理信息, 易于与其他地物类型和林地进行区分。 本研究采用灰度共生矩阵(gray-level co-occurrence matrix, GLCM)进行纹理特征构建, 根据植被光谱曲线, 植被在近红外B8波段计算纹理特征能够更好地区分植被等信息。 本研究将基于B8波段计算出17个纹理特征变量进行重要性分析, 为避免信息冗余, 只选取重要性排名前4的纹理特征变量。 由图4可选取均值和(sum of averages, SAVG)、 对比度(contrast, CONTRAST)、 差异方差(difference variance, DVAR)和相关性(correlation, CORR)4个纹理特征变量。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 基于B8波段的17个纹理特征变量重要性Fig.4 Importance of 17 texture feature variables based on B8 bands

地形特征反映了地表的形态和地貌特征, 提供了更加详细的地物环境背景与地物空间分布信息。 地形特征可以作为额外的输入特征, 与遥感影像数据一起用于训练分类器, 从而改善分类器的性能。 如表3所示, 本研究所构建的地形特征有坡向(Aspect)、 坡度(Slope)、 阴影分析(Hillshade)和高程(Elevation)。 Aspect 反映了地表在水平面上的朝向, 对地形起伏和地貌特征提供了重要信息。 Slope 表示地表在水平方向上的坡度程度, 对地形起伏和地势高低提供了直观的量化描述。 Hillshade 通过模拟光照效果, 呈现了地形表面的阴影分布, 增强了地形特征的可视化效果。 Elevation 是地表的海拔高度, 提供了地形的垂直信息。

表3

Table 3

表3(Table 3)

表3 特征变量 Table 3 Characteristic variable 卫星 特征类别 特征变量 Sentinel-1 SAR 极化特征 VV VH Sentinel-2A 光谱特征 增强型植被指数EVI 归一化差值建筑指数NDBBI 土壤水分指数SWI 大气阻抗植被指数ARVI 增强型建筑与裸地指数EBBI 水体提取指数WEI 纹理特征 均值和SAVG 对比度CONTRAST 差异方差DVAR 相关性CORR NASA SRTM 30m 地形特征 坡向Aspect 坡度Slope 阴影分析Hillshade 高程Elevation

表3 特征变量 Table 3 Characteristic variable

1.4.2 Sigmoid激活函数

由于不同光谱特征在特定条件下可能会产生非常大的值或非常小的值, 会造成未经归一化处理的光谱特征出现某些值域异常的特征值, 在训练过程中占据主导地位进而弱化其他正常特征值, 影响模型的训练效果。 在多光谱分析中, 不同指数的值域不同, 直接比较它们的值可能会导致误导性的结论。 本研究为解决此问题, 通过Sigmoid激活函数对研究中所用到的光谱特征进行归一化, 可以更容易地进行特征间的比较和分析, 提高模型精度。

Sigmoid激活函数适用于分类任务的输出层中, 输出值在0~1之间, 其还具有单调性、 平滑性以及中心对称性, 见式(1)

σ(x)=11+e-x(1)

式(1)中, x是输入值, 输入的值越大, 输出值也越大。 e是自然对数的底数(约为2.718)。 通过Sigmoid激活函数将研究中所用到的光谱特征值域映射到0到1之间, 将光谱特征标准化, 在后续分析和处理时更具有可比性。

1.4.3 时间窗口选取

由于不同植被物候特征、 生长强度以及对水分条件依赖的差异, 不同植被所表现出的光谱和极化的时序特征不同。 本研究采用Sentinel-1 SAR、 Sentinel-2A数据构建包括油茶林在内的五种土地类型的光谱和极化特征以选择最合适的时间窗口提取油茶林种植范围。

油茶对水的依赖性是识别油茶的关键特征。 因此, 研究中采用对水体敏感性较高且不受云雨天气影响的Sentinel-1 SAR VH极化特征, 用于检测土壤的湿度变化。 图5为五种土地类型在2021年的VH时间序列曲线, 图中油茶林VH时间序列曲线主要位于其他林地和不透水面之间且与不透水面有相交, 而在第23天和第316天之间, 油茶林VH时间序列曲线处于其他林地和不透水面之间且无相交。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 五种土地类型的VH时间序列曲线Fig.5 VH time series curves for five land types

鉴于油茶生长时土壤湿度较高的特点, EVI能够有效减少土壤背景噪声的影响, 提供更准确的植被健康和生物量信息。 研究中采用EVI作为监测油茶生长变化的依据, 构建油茶林2021年EVI时间序列曲线以区分油茶林与其他植被。 图6为五种土地类型在2021年的EVI时间序列曲线, 油茶林和其他林地的EVI时间序列曲线相似且与其他林地和耕地均有相交, 而在第22天和第141天之间, 油茶林EVI时间序列曲线处于其他林地和耕地之间且无相交。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 五种土地类型的EVI时间序列曲线Fig.6 EVI time series curves for five land types

根据五种土地类型所提取出的2021年五种典型地物的VH和EVI时间序列曲线, 选择时间窗口为2021年的第30天到2021年的第140天来研究吉安市油茶林种植范围。 其中在提取EVI时间序列曲线时采用Sentinel-2A的QA60波段对云和卷云覆盖进行剔除以减少异常值的干扰。

1.4.4 分类算法及精度评价

在机器学习中, 随机森林(random forest, RF)集成了多棵决策树, 使随机森林在分类中具有较强的抗过拟合能力, 能够有效地降低模型的方差。 梯度提升树(gradient boosted trees, GBT)通过组合多个决策树来构建一个强大的分类模型。 支持向量机(support vector machine, SVM)的关键思想是通过最大化间隔找到决策边界。 朴素贝叶斯(Naive Bayes, NB)是一种假设特征之间是相互独立的分类算法, 通过计算给定输入数据的每个类别概率, 并选择具有最高概率类别作为预测结果。

本研究所使用的交叉验证是评估机器学习模型性能的统计方法。 通过将训练数据分成若干个互不重叠的子集, 对每个子集轮流作为验证集, 其余部分作为训练集, 反复进行模型训练和验证, 获得对模型性能更准确的评估。 所有分类完成后, 采用四种混淆指标, 包括整体精度(overall accuracy, OA)、 Kappa系数、 用户精度(user accuracy, UA), 生产者精度(producer accuracy, PA), 用于分类结果的评价, 如式(2)— 式(5)。 整体精度是五种土地类型像元被正确分类到实际真实土地类型的概率; Kappa系数的作用是综合利用混淆矩阵中的参数进行离散性多元评价分类; 用户精度的作用是查看五种土地类型像元被正确划分到指定地类的概率; 生产者精度的作用是查看五种土地类型图斑与实际分类数据类型一致的概率。

OA(%)=∑i=1npiiN×100(2)

Kappa=N∑i=1npii-∑i=1n(pi+×p+i)N2-∑i=1n(pi+×p+i)(3)

UA(%)=piipi+×100(4)

PA(%)=piip+i×100(5)

1.4.5 卷积神经网络

(1) CNN网络

采用卷积神经网络(convolutional neural network, CNN)对油茶林进行分类。 设计并训练了一个深度CNN模型, 包括多个卷积层、 池化层和全连接层。 为保证深度学习与机器学习的数据源相同, 在研究过程中将分类所构建的所有特征均从GEE中导出后再导入深度学习模型中, 避免了因数据源不同而导致的结果差异。

(2) CNN-Bi-LSTM网络

为了进一步提高分类精度, 将卷积神经网络(convolutional neural network, CNN)与双向长短期记忆网络(Bi-LSTM)结合。 Bi-LSTM能够捕捉时间序列数据的前向和后向依赖关系, 将Bi-LSTM应用于时间序列遥感数据, 可以更好地理解和利用时空动态变化。 遥感数据具有时间跨度大、 数据量大的特点。 Bi-LSTM可以处理较长时间序列, 捕捉长时间段内的变化趋势, 增强模型的预测能力。 CNN与Bi-LSTM模块结合的模型图如图7所示。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 CNN与Bi-LSTM模块结构Fig.7 Structure of CNN and Bi-LSTM modules

1.4.6 余弦相似度

余弦相似度是一种衡量两个向量之间相似度的度量方法, 可以用于分析两个向量的方向是否相似, 如式(6)所示。 以2021年吉安市样本点的光谱向量为基础, 通过余弦相似度分析2021年吉安市样本点是否可以用于2019年、 2020年、 2022年、 2023年, 并将适合于其他年份的样本点保留并进行分类, 不适合的样本点则舍弃。

CosineSimilarity=A·B‖A‖‖B‖(6)

相似度阈值是一个预定义的数值, 用于判定两个样本点之间的光谱差异。 当计算的余弦相似度值大于或等于这个阈值时, 认为两个样本点的光谱不同, 即不是同一土地类型, 反之则认为它们相似。 通过对比2019年至2023年同期遥感影像, 共选取200个土地类型未变化的样本点, 每个土地类型含40个样本点, 并对每个样本点逐年采集研究中所涉及的各光谱波段, 得到最大和最小值范围^[16], 如表4所示。 结果表明, 油茶林的波动范围相对稳定, 各工地类型样本点在不同波段之间的波动范围处于0.035~0.534之间。 最终将相似度阈值的上限阈值设定为0.534。

表4

Table 4

表4(Table 4)

表4 无土地类型变化的样本点各波段阈值信息 Table 4 Threshold information of each band for sample points without land class change 波段 水体 不透水面 其他林地 耕地 油茶林 B2 0.101 0.509 0.036 0.041 0.040 B3 0.124 0.534 0.058 0.058 0.048 B4 0.159 0.528 0.035 0.089 0.085 B5 0.18 0.508 0.122 0.068 0.103 B6 0.228 0.459 0.290 0.179 0.104 B7 0.261 0.449 0.337 0.212 0.124 B8 0.223 0.493 0.348 0.293 0.131 B8A 0.268 0.429 0.352 0.233 0.133 B11 0.146 0.526 0.195 0.183 0.167 B12 0.198 0.514 0.102 0.164 0.145

表4 无土地类型变化的样本点各波段阈值信息 Table 4 Threshold information of each band for sample points without land class change

将分类场景分为两类, 一类为光谱特征经过归一化处理, 另外一类为光谱特征未经归一化处理, 两类分类场景均包含纹理特征、 地形特征和极化特征。 运用4种机器学习分类器对研究区进行油茶林识别, 得到研究区整体精度OA和Kappa系数以及油茶林用户精度UA和油茶林生产者精度PA。 由图8可知随机森林在4种机器学习分类器中对油茶林识别精度最高, 且在四种机器学习分类器中, 经过归一化后的光谱特征对油茶林识别精度明显高于未经归一化处理的光谱特征。 故本研究选用归一化光谱特征、 纹理特征、 地形特征、 极化特征和随机森林分类器作为油茶林识别最优方法来构建分类场景, 以进行后续的研究。

吉安市2021年油茶林面积分层图和空间分布如图9所示。 2021年吉安市油茶林总面积为184.488 1万亩, 其中被誉为“ 中国油茶之乡” 的遂川县油茶林种植面积最大, 占2021年吉安市油茶林总面积的27.67%。 种植分布呈现由地势较高向地势较低递减的趋势。 由图9可看出油茶林的种植普遍呈零散分布, 种植地点常见于家庭农户附近的山坡地和自留地, 利用土地的边缘区域进行栽培。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 不同机器学习分类器与不同分类场景的油茶林分类精度比较Fig.8 Comparison of the classification accuracy of Camellia oleifera forests with different machine learning classifiers and different classification scenarios

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 2021年吉安市油茶林面积分层、 空间分布图Fig.9 Stratified and spatial distribution of Camellia oleifera forest area in Ji'an, 2021

本实验深度学习框架为PyTorch2.2.2。 操作系统为Windows 11, CPU为AMD R5 5600 @4.0 GHz, GPU为NVIDIA RTX 4080super。 编程语言及版本为Python3.12.2, 学习率为10^-4, 迭代次数为200次, 批处理大小为64, 损失函数为交叉熵, 优化函数为Adam。

为保证模型之间的可对比性, 研究中CNN、 CNN和Bi-LSTM结合的深度学习模型采用与随机森林分类器一致的特征选取和样本数据, 分别对油茶林识别精度如图10所示。 CNN系列的深度学习模型对油茶林识别精度略低于随机森林分类器。 本研究只尝试了CNN系列模型, 并不能完全代表所有的深度学习模型对油茶林识别精度均低于随机森林分类器。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 不同模型对油茶林识别精度比较Fig.10 Comparison of identification accuracy of Camellia oleifera forests by different models

采用余弦相似度的2019年至2023年吉安市油茶林空间分布图如图11所示。 除2021年以外其他年份均以2021年样本点为基础所迁移得到。 由图11吉安市油茶林种植面积随年份的增加而呈现增加的趋势。

图11

Fig.11

	Figure Option ViewDownloadNew Window
	图11 2019年至2023年吉安市油茶林空间分布图Fig.11 Spatial distribution of Camellia oleifer forests in Ji'an, 2019— 2023

以2019年吉安市油茶林种植面积为基础的2020年至2023年吉安市油茶林种植面积增量如图12所示。 研究区内各县区的油茶林种植面积均在增加, 其中遂川县、 永丰县和永新县在2020年至2023年的整体增加速度较快, 吉州区、 青原区和吉安县在2020年至2023年的整体增加速度缓慢。

图12

Fig.12

	Figure Option ViewDownloadNew Window
	图12 2020年至2023年吉安市各县区油茶林增量Fig.12 Increase in Camellia oleifera forests in Ji'an counties and districts from 2020 to 2023