曹庆安,冷 鹏
(1.江西省核工业地质调查院,江西 南昌 330038;2.江西核工业测绘院集团有限公司,江西 南昌 330038)
水稻是我国种植面积最大、总产量最高的粮食作物,在粮食生产中具有举足轻重的地位[1]。水稻的年播种面积占粮食作物总播种面积的1/3[2],产区主要分布于长江中下游地区,耕作制度以一年两季为主[3]。目前,由于国内外环境的复杂多变,我国不断提高对粮食安全的重视程度,其中对于早稻播种面积的监测是预测全年粮食产量的重要环节之一,及时准确地监测早稻种植面积情况对全年粮食种植的宏观政策预判与调控具有至关重要的意义。
近年来,低空与高空遥感技术的发展突飞猛进,但卫星遥感技术因其数据获取范围广、时效性强、成本低等优势,广泛应用于农作物种植的监测[4-5]。目前,基于光学传感器获取的遥感数据已经难以满足农作物监测的需求,特别是南方种植区的多云多雨天气对光学遥感数据的获取造成了较大干扰。合成孔径雷达(Synthetic Aperture Radar,SAR)是一种有较强穿透性的高分辨率成像雷达,能够在复杂气象条件下进行全天候的对地监测,目前已广泛应用于自然资源、农林种植等领域[6]。
基于雷达的农作物种植监测研究已经成为农业资源遥感领域的热点方向[7],其中,水稻田因其具备高湿度的基底面,与旱生作物及其他地表覆盖类型具有较大的区分度,在特定物候周期内利用其后向散射系数的差异,可以实现较为准确的识别。然而,由于单时相数据的获取存在稳定性不强和对比度不高等缺点,研究者们尝试利用多时相、多极化的SAR数据提升识别精度[8]。目前,基于SAR的农作物识别方法大多是利用单一散射系数,但由于作物的下垫面散射与其生长周期存在较高的相关性,可能存在较大的误差并导致最终识别的精度不理想。
本文选取早稻不同生长期内连续时序双极化Sentinel-1A雷达数据,通过对比样本点间不同极化方式(VV、VH)后向散射系数的变化特征,研究不同影像分类结果的精度,初步提取作物的种植图斑。通过引入特定时间点的多光谱数据、坡度数据,对提取的作物图斑进行多因素决策分析,以期提升早稻种植图斑的识别精度,为早稻种植面积的监测提供快速、准确、便捷的方法。
高安市是江西省宜春市下辖县级市,位于江西省中部偏西北,地理坐标介于115°00′~115°34′E、28°02′~28°38′N,属长江中下游平原区,地形北高南低,为低丘岗地与河谷平原相间的地貌,属亚热带季风气候,年平均雨量1560 mm,降水期主要分布在4—7月份。高安市耕地面积约占全市土地总面积的30%,是江西省重要的水稻产区,种植制度以两季为主,部分地区也种植单季稻。
1.2.1 SAR数据 本文使用的基础雷达数据为2014年发射的Sentinel-1A卫星C波段的SAR数据[9]。该卫星的重访周期为12 d,主要包括条带(SM)、干涉宽幅(IW)、超宽幅(EW)和波(WM)模式4种数据。本文以干涉宽幅的双极化(VH、VV)数据为基础,时序范围覆盖了研究区早稻的整个生长周期。
1.2.2 地类样本数据 为精确获取研究区早稻的图斑位置和播种情况等信息,于2021年3月下旬—7月底针对不同的土地利用类型采集了典型的地类样本信息,包括早稻、蔬菜、林地、水面、建筑物以及其他地类等6类。此外,于2021年5月底—6月初通过无人机获取高分辨率航空正射影像,对研究区内的典型地类样本进行了二次筛选。最终确定了401个样本集,其中早稻132个、非粮农作物60个、旋耕田71个、林地30个、水面16个、建筑物60个以及其他地类32个。为验证像素分类的精度,采用分层随机选取了70%的样本作为训练样本,剩余30%的样本用于分类验证(表1),各类型样本点见图1。此外,2021年研究区的早稻物候期分别为:3月为播种期,4月为移栽(返青)期,5月为分蘖拔节期,6月为孕穗抽穗期,7月为灌浆成熟期,8月为收割期。
图1 研究区样本分布示意图
表1 地类样本统计情况
1.2.3 其他试验数据 本文用于分类决策的数据还包括:研究区30 m分辨率的ALOS WORLD 3D DSM数据,该数据由JAXA构建并向公众发布;多光谱数据,源自NASA与美国地质调查局(USGS)合作开发并发射的Landsat 8的星载OLI陆地成像仪,时序数据覆盖了早稻播种、生长和收割期;2021年度研究区的国土调查现状图斑数据。
1.3.1 SAR数据预处理 利用ENVI 5.6+SARscape软件对研究区的15幅SAR影像数据(2021年3月5日—8月20日)进行轨道数据精处理、辐射定标、多视处理、地形矫正以及滤波处理,经上述步骤后得到了多时相VV和VH时序影像的雷达后向散射系数[10]。
1.3.2 影像分类处理 本研究选取了3种常用的遥感监督分类方法对时序雷达后向散射影像图进行分类,并对分类结果进行精度评定,通过对比不同方法的总体精度与Kappa系数选取出最优的方法,最后对分类结果进行最小图斑处理并形成了初步的作物耕种图斑分布图。
1.3.3 土地利用图斑处理 本研究引入了最新的土地利用图斑,其能够反映研究区土地利用的类型,即实地地块的种植适宜度(是否为耕地)以及耕种属性(水田、旱地),可较为准确地提供适宜种植水稻的区域。因此,提取其中的水田图斑范围作为本次研究的早稻耕种本底数据能够极大地提高种植范围的提取精度。但土地利用现状图斑不能代表某一时刻的种植情况,还需要结合实地影像和其他参考资料进行综合分析。
1.3.4 植被指数 植被指数常用于监测植被的变化情况和绿色植物的分布情况,其能够有效反映研究区地表植被的生长状况。其中,归一化植被指数(NDVI)广泛应用于检测生物量、植被覆盖度、叶面积指数等植物生长状况指标[11-12],计算公式为:
式(1)中:Nir为近红外波段的反射值,R为红光波段的反射值。
NDVI能够直观地反映地表覆盖的量化特征。在早稻的生长期间,由于早稻的生长状况具有高度的一致性,其植被指数也具有极高的相似性。因此,通过对比样本点的NDVI差异,能够有效区分出早稻与其他植被覆盖样本。本研究利用早稻生长期内的播种初期、生长旺盛期2期以及收割期4期的NDVI灰度图,对相邻2期作差,获取了3幅NDVI差值灰度图。由于NDVI的系数一般介于-1~1,其差值可能存在负值。为突出显示NDVI的变化情况,需将作差后的结果统一归正后求倒数,该过程的计算公式为:
式(2)~式(3)中:NDVIΔ为相邻2期的植被指数的差值,NDVIi为第i期NDVI,NDVI′为植被指数差值归正后的倒数。
1.3.5 坡度数据 研究区位于长江中下游平原,主要是无坡耕地,基本分布在平原和河谷等地形相对平缓的区域。根据“三调”的地形坡度分布情况,耕地的坡度基本在0°~25°。通过30 m分辨率DEM获取研究区的坡度数据,计算公式为:
式(4)中,α为坡度,ΔH为高程差,D为水平距离。
1.4.1 最小距离分类法 最小距离分类法是利用训练样本对象对分类地物进行特征分析,计算出不同样本地类的均值和标准差向量[13]。一般是将待分类图像的像元与样本之间的均值向量的距离进行比较,数值最小的则归入该类样本。该分类方法中距离就是判别依据,距离越小说明距离该类别越近。距离的计算方法包括欧几里得度量和折线距离[14]。其中欧几里得度量的计算公式为:
式(5)中,N为波段数量,xi为第i个波段的像元值,Mij为第j类在第i个波段的均值。其分类的原则是将xi归入到Dej值最小的类别。
折线距离是计算像元值在各个波段中与类均值的距离差值的绝对值之和[15],计算公式为:
1.4.2 最大似然分类法 最大似然分类法又称为贝叶斯分类法,其面向待分类影像的波段数据的多维正态分布构造分类函数,首先统计每种样本地物的均值、方差和协方差等,然后将待分类像元代入构造好的分类函数中,计算其分类概率[16],计算公式为:
式(7)中:Si为第i类的第n个波段的协方差矩阵,i为波段数量,μi为均值向量,x为随机变量,P(Gi|x)为类别Gi对应x的概率。
1.4.3 SVM分类法 支持向量机(SVM)对小样本集具有良好的支持性,被广泛应用于遥感影像信息的分类提取[17]。SVM的原理是利用核函数将样本数据升维后重新规划,寻找最优惩罚因子和松弛变量。利用不同样本数据间超维度平面间隔大小作为分类依据。核函数一般使用径向基函数[18],公式为:
式(8)中:K(xi,xj)是SVM的核函数,γ为松弛向量。
影像分类结果的评价方法包含用户精度和制图精度[19]。精度评价存在一定程度的相对性,依赖于实际地类的结果验证又具有较大难度[20]。本研究利用遥感影像与其他多因素相结合的方法获取研究区的早稻种植图斑,并应用样本图斑对分类结果进行精度验证。
本文早稻种植图斑的提取精度是以分类结果落入样本地类的图斑面积占比进行评价。由于样本地类的面积大小不一致,其单个图斑的分类精度势必对总精度的影像存在差异,因此采用图斑面积权重的精度作为精度评价结果,计算公式为:
式(9)~式(10)中:E为精度,S为分类面积的总精度,Skn为单个图斑占总样本面积权重,Si为第i个样本图斑赋权后的精度,n为样本总数。
相对于监测其他类型的植物,农作物的生长周期具有固定性和可控性的特点。在特定区域内,由于水稻种植的时间也相对固定,该区域内水稻的生长状况也会具有高度一致性。因此,通过比较雷达数据的后向散射系数,可将早稻田与其他类型地物进行区分。
本研究的样本训练包括早稻、非粮农作物、旋耕田、林地、水面、建筑物以及其他地类等7个类别,按照早稻生长期的15个时相的VV和VH极化图像进行统计分析(图2)。由图2可知,不同类型样本在不同的极化图像和时间段内的后向散射特性具有明显的差异,VV极化下的散射系数比VH极化下的要高;在2种极化条件下,水体都呈现出最低的散射系数,而建筑物都呈现出最高的散射系数,但水体在VH极化条件下的散射系数最为稳定;林地、旋耕田以及其他非粮农作物的散射系数表现出稳定的波动性,且三者的趋势具有较高一致性,旋耕田由于地表裸露且含有大量的水分,其散射系数明显比其他典型地物低,且波动性小;VH和VV极化条件下的早稻样本的散射系数均呈现出了较为明显的波动性,如在4月下旬和7月下旬存在波谷,推测分别是处于播种期和收割期,5月中旬—7月上旬的散射系数呈现出明显的上升趋势,其反映出早稻处于生长的高峰期。从不同极化类型的地物分离度来看,VV极化不同地物散射系数区分度大,波动稳定性也相对较差,VH极化下不同地物间的散射系数差异较小,波动稳定性好。
图2 不同极化条件下的典型地物后向散射系数
以VV、VH极化为基础,根据早稻生长的不同周期,将波动的3个阶段分别以红、绿、蓝3种颜色显示,其中4月22日为红色、7月15日为绿色,8月8日为蓝色,形成了不同极化条件下的后向散射系数假彩色合成图像(图3)。
图3 VH极化(a)、VV极化(b)条件下的雷达后向散射系数RGB合成图
可分离度是衡量不同地物差异程度的量化指标,一般采用J-M距离(Jeffries-Matusita distance)对训练样本进行可分离性的定量评价[21]。J-M距离一般介于0~2.0之间,数值越大说明样本之间的可分离度越高,一般来说大于等于1.8可以较好地区分不同的样本。由表2可以看出,2种极化条件下的早稻样本数据与其他样本数据的分离度均超过了1.8,其中VH极化数据明显比VV的分离度要高。在VH极化条件下,早稻样本与非粮农作物的分离度相较于其他样本稍低,这是由于该时段一些水生作物或高湿度下垫面作物与早稻有较为相似的雷达反射特征。此外,4—7月份是全年降雨集中期,过高的水湿条件也对分离度的影响较大。因此,样本分离度不能将早稻与其他作物完全区分。
表2 不同极化条件下的时序数据样本可分离度评价表
采用最小距离法、最大似然法和SVM分类法对研究区VH极化条件下的3个典型时相地物后向散射RGB合成图进行分类(图4),并计算了分类后的总体精度和Kappa系数。由表3可知,SVM分类法的总体精度较高,达到了87.43%,Kappa系数为0.8240。
图4 VH极化条件下的分类结果图
表3 不同分类方法的精度评价
雷达的后向散射对早稻的提取可能存在一定的误差,本研究结合土地利用现状图斑、植被指数、地形坡度等多因素判定的方法对分类结果进行优化提取。
2.4.1 土地利用现状图斑 土地利用现状图斑是在经过实地调查的基础上形成的土地利用基底数据,其反映了该地块的土地利用类型,主要包括耕地、林地、园地、建设用地、未利用土地等,其中耕地继续细分为水田和旱地。本研究根据“三调”数据提取研究区的水田图斑(图5)。
图5 研究区土地利用现状水田图斑分布图
2.4.2 植被指数差值 对早稻生长过程的4期Landsat TM影像分别提取平均NDVI,各个影像的具体时间分别为2021年4月9日、5月27日和7月19日和8月20日,分别代表早稻种植过程的移栽(播种)期、抽穗期、灌浆成熟期和收割期。由表4可知,早稻在生长期的叶面生长状况存在明显的差异,表现为移栽(播种)期最低,灌浆成熟期和抽穗期最高;从样本平均NDVI差值来看,早稻图斑具有明显的波峰,并持续了2个监测时间点,直至8月20日才有明显下降;非粮农作物、其他地类以及旋耕田的指数差值变化不明显,且区分度不高;林地由于其植被生长特征,NDVI持续处于高值;建筑物、水面持续处于较低的水平。
表4 典型样本的平均NDVI
经过分析,拟将分类的训练样本地类归为:早稻、水面和建筑物、林地以及其他地类4类。为突出NDVI变化的情况,将依据4期影像的NDVI差值先正值化,然后分别求其倒数合成3幅灰度值图,分别赋红绿蓝合成假彩色图(图6)。采用SVM分类法对样本点进行分类,最终提取了早稻种植图斑。经精度评价,整体分类精度达到了89.64%,Kappa系数为0.8710,具有较强的可信度。
图6 NDVI差值时序影像合成图(a)及分类结果图(b)
2.4.3 坡度提取 通过DEM提取研究区坡度数据并对其进行分析。按照0°~5°、5°~15°、15°~25°、>25°等4个级别进行分级。最终形成了坡度分类结果(图7)。
图7 研究区地形坡度图
将后向散射系数灰度值分类结果图、土地利用现状图斑、植被指数差值分类结果图以及坡度提取结果图进行叠加分析,最终形成了研究区早稻种植分布图。根据图斑统计结果可知,2021年研究区早稻种植总面积为341.27 hm2,其主要分布于城西河谷平原地区(图8)。
图8 研究区早稻种植图斑分布
对最终提取的研究区早稻种植图斑与预留样本的图斑进行对比分析,精度评价将最终提取的早稻种植图斑落入预留样本图斑内的面积占比作为单个样本的分类精度,总体精度为所有预留样本的总面积的占比为权重计算单个样本赋权后精度进行汇总(表6~表7)。最终评定的分类面积精度为92.22%,赋权精度为94.76%。
表6 多因素决策的分类结果
表7 验证样本赋权分类结果
本文以Sentinel-1A的雷达影像数据为基础,对高安市的早稻种植面积进行了监测。以VH和VV极化方式对研究区早稻种植生长物候期的15期雷达影像进行了后向散射系数的分析,利用其变化情况组成时间序列影像。采用不同的分类方法对研究区的影像进行分类和评定,通过引入土地利用现状图斑、植被指数和坡度等多因素指标对分类结果进行优化,最终形成了早稻种植面积提取结果并进行了精度评价,得出以下结论:
(1)不同极化的雷达影像数据对不同地表覆盖类型的后向散射系数具有一定的差异。其中水体具有持续的较低的后向散射特征,建筑物和裸地表等具有持续的较强的后向散射特征。水稻的后向散射特征与其生长期基本匹配,主要表现出生长初期较低、生长期明显升高、收割期后显著回落的特征。林地以及其他植被覆盖具有较为稳定的后向散射系数。
(2)利用面向对象不同影像分类的结果发现,SVM的分类方法明显优于最大似然法和最小距离法。但最小距离分类方法对后向散射持续较高的建筑物和裸地等分类精度较优,最大似然法对后向散射较为稳定的林地等乔本植物覆盖分类精度较优。
(3)利用NDVI差值变化情况能够较好地将水稻等植被指数变化较为明显的农作物提取出来,本文提出了一种新的利用NDVI提取水稻图斑的方法:利用NDVI差值序列影像组合波段数据,基于样本地类对其进行SVM分类,最终提取的水稻种植图斑精度达到89.64%,能够较为准确地区分各个类型的地物。
(4)引入了土地利用现状图斑和地形坡度数据对提取结果进行优化,采用多因素分析的方式可以进一步提升水稻种植面积提取的精度。
(1)雷达后向散射系数在对农作物的识别中易受到基底水热条件的变化影响。一般来说,透水层反射的敏感度较大。本研究未能深入探讨不同地表覆盖在后向散射系数上的特征变化情况,这将是今后进一步研究的方向。
(2)旋耕田的雷达反射系数与早稻种植初期的和收割期后的水田反射特征具有较高的相似度。此外,其他水生作物(如茭白、莲藕)等也会对早稻种植面积的提取造成一定的影响。虽然本研究为排除旋耕田和水生作物对分类结果的负面影响进行了一定的探索,但并未深入进行定量统计和实地验证,今后需要在该方向上做更多的探索。
(3)基于landsat TM影像的植被指数虽然能够反映不同植被的生长状况,但其分辨率较低(30 m),采用更高分辨率的影像如Sentinel-2、GF-2等能够提供更优的数据支撑。
(4)利用本方法对国产的GF-3卫星雷达影像进行早稻种植面积提取可以作为本方法的衍生研究。