1. 引言
近年来，随着我国农业结构的改变与城市化进程的加快，使得部分地区出现耕地“非粮化”倾向，致使部分耕地被违规改种草皮等非粮作物，对粮食安全造成威胁[1]。草皮因其生长周期短、经济效益高，种植的范围越来越广泛(如图1所示)，中国草皮产业规模已列世界第三位，仅次于美国和欧盟[2]。但是，占用耕地种植草皮的传统生产方式对耕地土壤危害严重，平均每次种植会带走2.5~5厘米的表层土壤，会逐步带走表层土壤，破坏耕作层；造成养分流失，土壤有机质得不到补充，导致土壤贫瘠、沙化[3]，因此对耕地内种植草皮进行大范围、快速、精准的监测十分关键。
卫星遥感已经是当前我国开展大范围耕地种植情况监测的主要手段之一，但当前却少有使用其进行识别耕地内种植草皮的研究。合成孔径雷达(SAR)影像技术因其能够穿透云层、不受天气条件限制，成为监测耕地内种植情况的理想选择，通过SAR影像的后向散射变化，可以获取耕地内作物种植的动态变化信息，从而实现对耕地使用情况的快速监测。本文将探讨基于时序SAR影像的岭南地区耕地内种植草皮的识别方法，旨在为耕地保护和合理利用提供科学依据和技术支持。
Figure 1. Photo of turf grass on arable land
图1. 耕地上种植草皮的照片
2. 研究现状
使用遥感进行耕地内草皮提取的研究很少，陈艺华等利用高分多光谱影像、使用基于田块对象提取了田块纹理特征进行了镇级的全域草皮提取，取得了较好的效果[4]；席瑞等使用了多时相植被指数和基于对象的方法对稻田内草皮进行了提取，主要利用了水稻的植被指数与草皮植被指数动态变化特征加以区分。这些研究方法均依赖于多光谱影像的获取频次和成像质量。
但要应用于南方地区尤其是岭南地区，多光谱影像主要存在两个局限性：一是易受云雾遮挡影响。由于我国多雨多云地区的雨热同期的特性，草皮种植期内往往大气窗口狭窄，难以获取到无云污染的高分辨率多光谱影像；二是耕地内草皮在可见光影像上易和一些更常见的地物容易混淆。耕地中的生长草皮往往呈现较均匀的绿色，容易与水稻、玉米、自然草地等混淆，而收割之后的草皮呈现土壤的棕黄色，容易与裸土、翻土后的耕地混淆。
近年来，随着越来越多大范围、高频次的合成孔径雷达(Synthetic Aperture Radar, SAR)数据可用，研究者逐渐转而使用SAR进行农作物监测[5]。SAR通过宽带信号脉冲压缩以及虚拟孔径合成来获得高分辨率图像，且不受时间和天气的影响，全天候全天时工作，同时其获取影像的纹理、结构信息丰富，可视为光学遥感数据的有效补充。并且SAR系统支持多种极化方式(如HH、VV、HV、VH)，可以提供额外的地表信息，有助于提高对地物的识别能力。
但目前尚未有文献研究利用SAR对耕地内种植草皮进行监测，因为单时相SAR数据在耕地识别应用中也有其局限性：一是其在使用时仅依据后向散射信息或者极化信息，对于地物的区分度有限；二是田块区域地物破碎化情况较突出，田块分布密集，导致耕地作物与草皮的边界识别不够清晰。对于面积较小、形状不规则的田块，覆盖范围广但分辨率相对(多光谱)较低的SAR影像混合像元问题严重，叠加上SAR本身难以去除的乘性噪声，使得其较难满足田块级别的精细监测需求。
由于近年来短重返周期遥感周期的影像(如哨兵系列卫星)更易获取，使用时序SAR进行的作物识别方法十分流行，尤其是水稻等不同物候期后向散射变化明显的作物，不同作物生长周期内，耕地内水分变化、植株的生长为不同种类的作物刻画了不同的时序曲线[6]，这成为区分不同地物的重要依据。基于此，本论文将探索基于时序SAR影像的草皮识别方法的可行性。
遥感影像地物分类与识别常用的算法包括机器学习算法和深度学习算法[7]。深度学习方法[8] [9]能够从卫星图像数据中检索复杂模式和信息特征，能够获得很好的性能。然而，基于深度学习计算机视觉的方法多依赖于可见光数据，且对于目标正样本的草皮样本数量要求较高，同时，基于深度学习的方法计算成本较大，推广模型至较大的尺度时对GPU等硬件要求较高。最近的研究表明，机器学习可以用较少训练数据集处理学习任务，与深度学习相比，其结果却不相上下，支持向量机SVM [10]-[12]、随机森林RF [13]和XGBoost [14]等机器学习方法因其计算复杂度较低、可解释性相对较高而受到研究人员的青睐[15]。其中SVM在处理高维、噪声数据集时容易出现过拟合，同时忽略了筛选和评估特征，其二分类特性也限制了其在遥感领域的应用[16]；RF在某些噪音较大的分类或回归问题上会存在过拟合问题，同时对于有不同取值的属性的数据，取值划分较多的属性会对随机森林产生更大的影响。而XGBoost算法由于具备运算高效、灵活性强、预测精准、便于并行计算等特点，近些年在地理信息挖掘、图像模式识别领域得到广泛应用[17] [18]。基于此，本论文将探究使用XGBoost算法进行基于时序SAR影像的草皮识别。
3. 研究方法
3.1. 田块分割
基于特征的方法一般有两种常用预测单元，基于像素和基于对象。基于对象的遥感影像分类方法通过将图像分割成多个相对独立的地物对象，考虑对象的形状、纹理和空间关系等特征，从而提高了分类的准确性和鲁棒性。相较于基于像素的方法，基于对象的方法能够有效减少噪声干扰，并更好地捕捉地物之间的上下文信息。
耕地的田块通过Xie等[19]提出的方向引导后处理的耕地边缘检测方法，该方法首先利用HRNet和OCR模块进行边缘检测，以获取更精确的农田边缘信息；其次，引入ConAttn模块生成方向信息，并通过断点连接算法连接断开的边界线，从而得到完整闭合的农田边缘；最后，将闭合的边缘线转换为多边形结果，生成精细的农田地块。该方法有效解决了传统语义分割方法难以提取精细农田地块的问题，并能够生成更完整、更详细的农田地块。
基于此方法提取的田块，以田埂为单元边界，比地块单元更精细，以稍高的计算代价，解决了地块破碎地区的地块内种植属性异质性高(同一地块多块田不同种植属性)的问题。
3.2. 基于田块对象的特征提取
在获取到研究区耕地内所有田块后，可以对田块对象内的像元进行的特征提取，主要包括像元直接统计特征、像元直方图统计特征和纹理特征，不同影像所用的各统计特征如表1。
Table 1. Parcel based image features
表1. 以田块对象为单元的影像特征
影像类型 |
像元统计 |
直方图统计 |
纹理特征 |
可见光 |
最大值、最小值、均值、众数、中位数、像元个数、离差、第10百分位数、第90百分位数、标准差 |
偏度、峰度 |
灰度共生矩阵同质性、能量、相关性 |
SAR |
最大值、最小值、均值 |
|
|
统计特征指像元的直接统计值和像元直方图的统计值。
像元统计值是指以田块为单元的所有像元的统计值，包括最大值、最小值、均值、众数、中位数、像元个数、离差、第10百分位数、第90百分位数、标准差。其中最大值、最小值、均值分别表示田块内像元的光谱反射率的最高、最低值，以及平均反射率。这三项特征用于反映地物的整体反射强度；众数、中位数、像元个数用于描述田块内像元反射值的集中趋势；离差度量田块内各像元的光谱反射值与均值的偏离程度，表明地物表面反射率的均匀性；第10百分位数、第90百分位数相比于最小值和最大值更能对抗极端的异常值，可以减少噪声的影响；标准差衡量田块内光谱像元值的离散程度。像元直方图统计值是对像元数据的更高阶矩的表征，偏度描述像元值分布的对称性，峰度描述数据分布的尾部厚度和峰的尖锐程度。
纹理特征是用来描述图像中像素灰度值在空间上的分布规律和结构特性的，其能反映图像的表面特性和视觉模式。灰度共生矩阵是常见的描述纹理的统计方法[20]，其统计值包括灰度共生矩阵的熵、同质性、能量、相关性、ASM、不相似度等。因为RGB各波段纹理特征几乎相同，为了减少特征冗余，灰度共生矩阵的选择了RGB合成(0.299 × Red + 0.587 × Green + 0.114 × Blue)的灰度波段和近红外波段。灰度共生矩阵计算时选用了2个方向(90˚和180˚)和3个距离(1、3、5个像元)。由于灰度共生矩阵的特征部分统计指标之间也有较大的相关性，故灰度共生矩阵统计指标仅选择同质性、能量、相关性三个指标，以减少特征的冗余性。
因为研究区耕地田块破碎，部分小田块面积较小(小于400 m2)，对于分辨率为亚米-米级的可见光影像来说，上述统计特征和纹理特征均有意义。而对于分辨率为10 m的哨兵1-SAR影像来说，标准差等复杂的统计值和纹理特征失去意义或无法计算，故本研究仅使用基本像元统计值，即最大值、最小值和均值。
为了验证时序SAR影像特征的有效性，这里选用了一期光学影像进行包括像元统计值、像元直方图统计值和灰度共生矩阵统计值等三部分的光学影像特征(后简称为光学影像特征)的提取以进行控制实验。
3.3. 基于XGBoost的草皮识别模型
XGBoost (Extreme Gradient Boosting)是一种高效的梯度提升算法[21]，广泛应用于分类和回归问题。它通过组合多个弱学习器(决策树)形成一个强学习器，每次迭代根据前一轮的误差来调整模型。XGBoost通过使用正则化来防止过拟合，并通过并行处理和近似算法提高计算效率，适用于大规模数据集。XGBoost能处理缺失值能力：能够自动处理缺失数据，而不需要特别的数据预处理，这给实际遥感影像进行监测时偶有缺测的问题提供了一个简单快捷的解决途径。
用XGBoost算法对样本统计属性特征直接进行样本分类，为了减少模型过拟合情况，本研究不对其进行复杂的超参数网格搜索，超参数只调节了树的个数，选取范围是10、20、50、100、300、500、1000，其它超参数均使用XGBoost Python库(v1.7.6)的默认值。
4. 实验结果与分析
4.1. 研究区与数据
4.1.1. 研究区概况
本研究从岭南地区选取了广东省的两个区县级区揭西县和新会区进行实验，如图2。揭西县隶属于广东省揭阳市，位于广东省东部，山地和丘陵居多，平原仅占约14%，属亚热带季风气候，气候湿润，温暖多雨，光照充足，年降雨量2000毫米左右。新会区隶属广东省江门市，位于珠江三角洲西南部，西江和潭江下游，平原面积约占44%，属亚热带海洋性季风气候，雨量充沛，日照充足，常年温和湿润，无霜期长，年平均降雨量1800毫米左右。
两个研究区都是典型热带亚热带季风气候，水热充足，植物生长旺盛，尤其5~9月多云多雨，较难获取高分辨率可见光影像。耕地内均以水稻为主要粮食作物，容易与耕地种植的草皮产生混淆。另外，研究区丘陵山地较多的地方，地块破碎，种植复杂，往往一年种植多造作物，给耕地内种植物的识别带来了更大的挑战。
Figure 2. Study area of Jiexi County and Xinhui District
图2. 研究区域揭西县和新会区
4.1.2. 实验数据
鉴于研究区内多云多雨的特性，本研究中选取了可以可全天候观测的SAR影像作为识别的主要影像。研究分别获取揭西县、新会区两个研究区2022年1~10月的多时相哨兵-1A星1月~10月的C波段Level 1干涉宽幅(IW)观测的GRD数据，该数据重返周期为12天，共25景，空间分辨率为10米，经过ESA SNAP软件进行了轨道文件校正、热噪声去除、辐射定标、斑点噪声滤波(Lee-Sigma)和地形校正，得到研究所需的VH后向散射系数(dB)数据。
同时，为了进行对比，也选取了国产高分光学影像进行了测试。分别获取揭西县、新会区两个研究区同期的一幅干净无云高分辨率遥感影像，分别为2022年9月的高分六号和2022年3月的高分1B四波段多光谱(RGB+近红外)高分影像，分辨率为3米和2米，除了把光学影像用于特征提取之外，研究区的光学影像还用于田块提取程序，以提取研究区田块单元。
本研究中，种植属性分类样本数据从实地田野调查采集获得，以田块为单元，揭西县和新会区分别获取了11,208和3815个田块样本，分类占比如图3所示，为保证样本的完整性和代表性选取其中400平方米以上的图斑(即至少4个哨兵-1A SAR影像像元)，最终得到揭西县3747个和新会区3311个田块样本数据。在机器学习分类中，田块样本数据以7:3分层抽样划分为训练集和验证集。
Figure 3. Sample distribution map of Jiexi County and Xinhui District
图3. 揭西县和新会区田块样本分布图
4.2. 草皮时序特征曲线
在水稻分类的研究中，常用时序SAR的曲线来刻画其种植和物候特征[6]，为了研究草皮与其的时序特征区分，本研究绘制了耕地内三种地类的时序曲线，如图4为揭西县2022年耕地范围内田块样本SAR影像均值特征的分类平均及其0.5倍标准差时序曲线。
Figure 4. Classification average (solid line) and 0.5 times standard deviation (shaded area) time series plot of VH mean for field samples within the arable land area
图4. 耕地范围内田块样本VH均值特征的分类平均(实线)及其0.5倍标准差(填色)时序图
从均值水平上来看，水稻的后向散射系数均值较高，尤其在生长阶段和成熟阶段，作物表面的粗糙度和密度增加，导致电磁波反射增强，后向散射系数升高；而草皮的均值始终较低，大致与水稻的灌水时期相当，这可能和草皮的种植时的频繁人工浇灌相关。
从波动幅度上来看，代表水稻的绿色曲线随时间波动明显，反映了其在生长过程中的物候变化，随着植株的生长、成熟、收割，株高、叶片数量和含水量等因素均会改变作物的表面反射特性，最明显的特征是两个波谷，代表了双季稻的两次灌水期；而代表草皮的蓝色曲线变化波动幅度相对较小，但波动频率相对较高，这与耕地内种植草皮生长迅速、一年多次种植的特性相符。
通过对耕地范围内不同地物类型的SAR影像时序曲线分析，可以清晰地发现水稻、草皮和其他地物在后向散射系数曲线上的差异。时序特征曲线的差异是机器学习算法能够区分草皮与水稻及其他地物可分性的基础。
4.3. 分类结果及精度评价
两个研究区使用单期光学影像特征、时序SAR影像特征及其组合的分类实验结果如表2、表3所示。为了综合评估模型性能，使用总体精度和Kappa系数衡量模型在多类别中的分类性能，使用准确率、查全率和F1得分衡量草皮单类的分类性能。
Table 2. Performance of turf recognition in Xinhui District
表2. 新会区草皮识别精度表
组合类型 |
最优模型 |
总体精度 |
Kappa系数 |
准确率 |
查全率 |
F1得分 |
单期光学影像 |
300树 |
86.4% |
0.746 |
85.9% |
52.4% |
0.651 |
时序SAR影像 |
50树 |
94.2% |
0.893 |
94.7% |
84.8% |
0.894 |
单期光学影像 + 时序SAR |
50树 |
95.3% |
0.914 |
96.7% |
84.8% |
0.904 |
Table 3. Performance of turf recognition in Jiexi County
表3. 揭西县草皮识别精度表
组合类型 |
最优模型 |
总体精度 |
Kappa系数 |
准确率 |
查全率 |
F1得分 |
单期光学影像 |
1000树 |
79.9% |
0.624 |
81.7% |
46.2% |
0.59 |
时序SAR影像 |
50树 |
88.4% |
0.787 |
90.3% |
79.2% |
0.844 |
单期光学影像 + 时序SAR影像 |
100树 |
82.5% |
0.946 |
96.7% |
82.1% |
0.879 |
结果表明，仅使用时序SAR影像特征的XGBoost模型在两个地区的多类的性能以及草皮单类的分类性能上均展现出了明显优于单期可见光影像特征的结果，其中F1得分在新会区和揭西县分别为0.894和0.844，远高于仅使用单期光学影像特征的0.651和0.590。经田间调查发现，在耕地内种植草皮一般是散户耕种模式，研究区内相近的草皮种植区域在调查时可能是生长期或者刚收割结束，这导致其在光学影像上呈现出绿色和棕黄色等迥异的两种状态，这两种状态分别容易与生长期的水稻和裸土误分，所以识别的结果依赖于光学影像拍摄的时段是否刚好捕捉到关键的生长期或收割期。而本研究中，时序SAR数据共有25期数据，时间分布较为均匀，以完整作物生长时间序列曲线，时序SAR在连续观测和对地表水的敏感性使其在草皮与水稻、耕地内其他覆盖的区分上有较大优势，因而具有更好的表现。
同时使用单期光学影像特征和时序SAR特征的组合策略，均提升了识别精度，两个研究区的F1得分分别达到了0.904和0.879，但是提升幅度均不大，新会区提升了0.1，相比之下基数较小的揭西县提升达到0.35。单期光学影像特征相较于时序SAR特征，具有更丰富的统计特征和纹理特征，二者结合增加了各地物间的光谱差异性和可分性，这使得草皮与其他地物之间的差异更明显，识别精度更高。
综上可以看出，单独使用SAR影像已经能够达到较高的识别精度，足以满足草皮识别的实际需求。且在实际应用模型时，岭南地区的5~9月间无云的光学影像可遇不可求。
使用时序SAR影像数据统计属性特征训练好的模型对揭西县和新会区草皮进行识别，获取最终草皮的识别结果，草皮预测结果如图5所示。
Figure 5. Turf extraction result samples (Red for turf grass; Green for paddy rice)
图5. 草皮提取结果示例(红框：草皮，绿框：水稻)
5. 结论与展望
本研究针对岭南地区多云多雨及耕地破碎化等条件，提出了一种基于时序SAR影像和XGBoost算法的草皮识别方法，该方法利用基于田块的时序SAR，构建了推理速度极快的机器学习模型，同时实现了较高的识别精度。结果表明：(1) 草皮后向散射时序曲线与水稻等耕地内其他地物具有明显区别，其VH时序均值低于水稻，振荡频率高于水稻，体现了其种植特性；(2) 基于时序SAR特征的XGBoost模型验证结果表明，在新会区和揭西县的独立试验查全率和准确率均在80%和93%以上，表明模型具有较好的鲁棒性和可推广潜力；(3) 不同特征组合下的控制试验表明，仅使用时序SAR特征可以较准确区分光学上易混淆的草皮、其他草地、水稻等地物，而单期光学和时序SAR组合的模型取得最高精度。
本研究构建的XGBoost模型具有运行快速、抗数据缺失能力强的特点，易向多云多雨的华南地区推广应用，研究结果可为耕地“非粮化”精细监测和科学管控、保障粮食安全提供数据支撑和决策支持。
尽管本研究取得了较好的结果，但仍存在一些不足之处。首先，本研究主要依赖基础的哨兵-1A SAR影像的VH极化后向散射的基础统计特征，VV极化和极化SAR的特征如极化特征和相干性特征的应用尚有待进一步探索。其次，本实验的基础分类单元是基于光学影像分割的田块单元，结果受到田块分割算法效果的影响，如果耕地内田块漏分，则会导致相应的草皮漏分。
随着更高时空分辨率SAR数据普及率越来越高，如陆探一号、COSMO-SkyMed等，这些数据的引入是否能进一步提高耕地内种植草皮的检测准确性，这是下一步可以研究的方向。
基金项目
广东省科技计划项目(2021B1111610001, 2021B1212100003)，广东省自然资源科技项目(GDZRZYKJ2023001, GDZRZYKJ2024002)。
NOTES
*通讯作者。