基于时间序列全极化合成孔径雷达的水稻物候期反演

李宏宇，李坤，杨知

（1.上海市测绘院，上海 200063；2.中国科学院空天信息创新研究院，北京 100101；3.中科卫星应用德清研究院/浙江省微波目标特性测量与遥感重点实验室，浙江 湖州 313200；4.中国科学院大学，北京 100049；5.中国电力科学研究院有限公司，北京 100192）

水稻作为人类主要食物来源之一，世界上超过一半的人口以水稻为主食。我国是水稻种植和消费大国，水稻总产量约占全国粮食总产量的33.75%[1]。水稻在经济发展和社会稳定等方面具有举足轻重的地位。

遥感技术因其能大面积同步观测、重访周期短、时效性高等优势渐渐取代原始的人工测量方法，成为水稻生长监测的主要技术手段之一。目前，基于光学数据和全极化合成孔径雷达（synthetic aperture radar, SAR）数据对水稻长势监测、估产和生长参数反演等方面的研究已趋于成熟[2-6]。然而，由于水稻长时间生长在云雨天气下，云雨等因素会严重影响光学遥感数据质量。因此，利用光学遥感数据难以实现对水稻关键物候信息的准确提取。而合成孔径雷达因其全天时、全天候的特点，克服了光学数据易受天气影响的不足。早期由于雷达技术发展水平有限，因而多使用单极化或多极化数据进行水稻识别及长势监测研究[7-10]。随着雷达遥感数据发展至全极化时期，全极化SAR数据在水稻长势监测中的应用变得更为广泛。相较于单、双极化数据，全极化SAR 数据具有地物完整的极化矩阵、几何结构细节等信息。大量研究表明，全极化SAR 数据在水稻识别和物候期提取等方面更具优势。LOPEZ-SANCHEZ 等[11]利用TerraSAR-X 双极化数据进行水稻物候期识别研究，提取出了幼苗早期、分蘖末期、拔节期等5 个水稻物候期，识别精度在80%以上；PACHECO 等[12]使用TerraSAR-X 和RADARSAT-2 数据进行水稻物候期识别且取得了更高的精度。然而，这些研究对全极化SAR数据仅局限于使用数据强度特征、波段组合特征，并未充分利用极化信息，且由于全极化SAR数据中包含众多极化特征参数，鲜有研究构造对水稻物候变化敏感的有效表征指数。为此，本文基于多时相全极化RADARSAT-2数据，对研究区内水稻种类籼稻和粳稻进行识别；在此基础上，筛选出对水稻生长物候变化敏感的极化特征参数，通过分析水稻生长过程中生理结构变化导致的雷达响应特征，构造了一个能反映水稻生长变化特征的雷达物候指数（radar phenology index,RPI），分析水稻生长过程中时间序列RPI 变化特征，利用导数法提取时序RPI 曲线关键特征点（拐点、最值点），从而提取出水稻生长关键物候期。

1 研究区和数据源

1.1 研究区概况

本文以江苏省淮安市金湖县周边为研究区（图1），其是中国主要水稻种植区之一[13]。研究区的整体占地面积约为40 km×30 km，中心点经纬度坐标为33°07′05″N，118°59′55.14″E。研究区内水稻的种植模式主要为一年一熟；水稻种类主要分为籼稻和粳稻2种。籼稻的全生育期约为120 d，粳稻的全生育期约为150 d。研究区内除了水稻之外，还包含其他5种典型地物：森林、水体、城市、蟹塘和浅滩。

1.2 全极化SAR 数据源与地面同步测量实验数据获取

1.2.1 多时相全极化SAR 数据

本文获取了2015 年共6 景RADARSAT-2 全极化SAR数据，具体参数见表1。RADARSAT-2卫星的工作频段为C波段，其极化方式为全极化，重访周期为24 d。本文获取的SAR 数据均匀分布在水稻的全生育期内（播种、出芽、幼苗、分蘖、拔节、孕穗、抽穗扬花、乳熟、蜡熟和完熟阶段）。

1.2.2 地面同步测量实验

为了确保获取到的全极化SAR 数据能准确地反映出研究区内植被的生长状况以及其他地面信息，需要在卫星过境的同时进行野外地面同步测量实验。

图1 2015年9月16日RADARSAT-2全极化SAR数据RGB假彩色合成图Fig.1 RGB false color composite image of RADARSAT-2 full-polarization SAR data on September 16,2015

表1 获取的RADARSAT-2全极化SAR数据参数Table 1 Acquired RADARSAT-2 full-polarization SAR data parameters

在研究区内分别选择12块籼稻和19块粳稻样本田，每块样本田的面积都在120 m×120 m以上，能够保证每块样本田包含足够多的像元来进行极化特征参数的提取及分析。水稻地面同步测量实验的获取参数包括：水稻样本田的经纬度坐标、区域矢量文件及水稻物候信息等。此外，还对研究区内森林、水体、城市、蟹塘和浅滩5种典型地物类型进行调查取样，获取各样本的中心经纬度坐标和矢量文件。

1.3 多时相全极化SAR 数据处理

1.3.1 多时相全极化SAR 数据预处理

全极化SAR 数据预处理使用欧洲航天局（European Space Agency, ESA）开发的SNAP 5.0 软件完成。SAR数据预处理主要包括辐射定标、几何校正和相干斑滤波处理。

首先进行SAR 数据辐射定标处理。将SAR 图像的灰度值相对于标准雷达截面积进行定标，转化为后向散射系数。然后进行图像几何校正处理。雷达图像的成像方式为侧视成像，获得的图像为斜距图像，存在多种类型的几何畸变以及叠掩、阴影、透视收缩等现象。本研究中几何校正采用SRTM-90m 的数字高程模型数据，像元重采样的大小为5 m×5 m；在前2 步完成的基础上最后进行滤波处理。原始SAR图像中存在大量相干斑噪声，使得图像数值不能真实地反映地物的散射特性。因此，进行降噪处理是定量反演前的必要过程。本文采用增强型Lee 滤波方法去除SAR 图像中的相干斑噪声，滤波窗口大小为5×5。

1.3.2 多时相全极化SAR 特征参数提取

为了探究研究区内典型地物极化散射特性随时间的变化特征，利用野外实验获取的粳稻、森林等地物的样本矢量文件，基于2015 年多时相RADARSAT-2 数据进行全极化SAR 特征参数提取。所用到的软件为欧洲航天局开发的PolSARpro 5.1.2软件和SNAP 5.0软件，共提取了119个极化特征参数，其中部分参数如表2所示。

表2 全极化SAR部分特征参数提取列表Table 2 Partial characteristic parameter extraction list of fullpolarization SAR

2 基于时间序列全极化SAR 的水稻种类识别

一个完整的水稻生育期主要包括播种、出芽、幼苗、分蘖、拔节、孕穗、抽穗扬花、乳熟、蜡熟和完熟阶段。2015 年SAR 数据的获取日期对应的水稻物候情况如表1所示。

将所有地物样本中的2/3 作为训练样本，1/3 作为验证样本。利用训练样本矢量文件提取出多时相极化特征参数，统计训练样本极化特征参数的平均值并作为该地物在该时相下的极化特征参数数值。图2 表示研究区内典型地物2015 年部分极化特征参数时序变化情况。利用不同时相下典型地物的An&Yang 体散射分量（An&Yang3_Vol）、反熵（anisotropy）、水平极化与垂直极化分量比值（δHH/δVV）、香农熵（Shannon entropy）等极化特征参数数值之间的差异，将籼稻和粳稻从其他地物中识别出来。

香农熵可以定量地反映遥感信息的分散程度。粳稻和籼稻全生育期长短不同且生长起始时间相近，当粳稻处于抽穗扬花期时，籼稻处于乳熟期。水稻处于抽穗扬花末期时，体散射绝对占优，二次散射和面散射分量较小，使得香农熵处于较低的水平；水稻处于乳熟期时，水稻冠层不再封闭，对微波能量的衰减作用稍有减弱，导致体散射、面散射和二次散射分量均稍有上升，香农熵数值升高。因此，可以利用这一时期的香农熵将籼稻和粳稻区分开（图2）。

本研究选取决策树分类器进行水稻种类识别。图3所示为本文所使用的决策树算法。最终得到的识别结果如图4 所示。从中可以看出，粳稻主要分布在研究区的北部和东部，籼稻主要分布在南部和东部。

根据最终的分类结果以及验证样本得到混淆矩阵（表3），并进行分类精度评价。从中可以看出：所有地物的总体分类精度为94.73%，Kappa 系数为0.93。水体、蟹塘、浅滩、城市、森林、籼稻和粳稻的用户精度均在85%以上，制图精度均在80%以上，得到了精度较高的分类结果。

3 基于时间序列全极化SAR 的水稻物候期监测

3.1 水稻时序雷达物候指数曲线提取

3.1.1 雷达物候指数构建

图2 典型地物极化特征参数变化特征Fig.2 Change features of polarization characteristic parameters of typical ground objects

图3 水稻种类识别的决策树算法Fig.3 Decision tree algorithm for rice species identification

图4 研究区地物识别结果图Fig.4 Result map of ground object recognition in the study area

表3 基于多时相全极化SAR数据的研究区地物分类混淆矩阵Table 3 Confusion matrix of ground object classification in the study area based on multi-temporal full-polarization SAR data

选取Cloude-Pottier分解中的λ1分量、Yamaguchi四分量分解中的体散射（Yamaguchi_Vol）和面散射（Yamaguchi_Odd）分量来构造一种能反映水稻生长物候变化特征的RPI，其表达式见式（1）。利用体散射与面散射分量进行比值处理。由水稻的雷达响应机制可知，水稻的面散射与体散射之间呈近似反比的关系，即当水稻的面散射贡献较多时，体散射分量相对较小，其比值可以突出2 种散射机制之间的差异，结合λ1的强度信息，可以反映水稻在整个生长过程中主要强度的变化情况。因此，可以结合RPI 的时间序列数据特点识别出水稻关键物候期。本研究基于2015年6景SAR数据，选取19个水稻样本田进行RPI提取。

3.1.2 基于Savitzky-Golay（S-G）滤波的RPI 时间变化曲线重构

3.1.2.1 数据插值处理

由于RADARSAT-2 卫星数据的重访周期为24 d，且水稻的生长速度相对较快，导致每景卫星数据间隔时间内会错过某些水稻关键物候期。因此，本研究选择3次样条插值法对2015年水稻时序RPI数据进行插值处理。经3次样条插值后的数据曲线更为光滑、稳定，可以确保时序数据二阶导数的连续性。本研究经3次样条插值处理后的数据时间间隔为12 d。以第10 号样本田为例，插值前后的RPI时序数据对比如图5 所示，插值后的RPI 曲线取得了一定程度的平滑效果，保持了相对连续性，拟合出了未获取到卫星数据时相下的RPI，有利于后续准确反演水稻关键物候期。

图5 第10号样本田2015年水稻RPI插值处理前后对比图Fig.5 Comparison chart before and after RPI interpolation of rice in sample field No.10 in 2015

3.1.2.2 S-G滤波重构

时间序列RPI 仍会受到环境等因素的影响，存在数值过高或过低的现象，因而需要对时间序列RPI数据进行滤波降噪处理。S-G滤波法是一种基于最小二乘的卷积拟合算法，由SAVITZKY 等[28]在1964年提出，其在对有关植被参数进行滤波重构时不受时间和传感器等因素的限制[29-30]。本研究经过多次实验，设定拟合多项式阶次为2 和平滑窗口的宽度为5。重构后的时序RPI 如图6 所示。通过对比可以看出，经过S-G 滤波后的时序RPI 在没有改变基本形态的前提下更加平滑，部分区域的噪声被修正。

3.2 基于时序RPI 变化曲线的水稻关键物候期识别

3.2.1 水稻时序RPI 时间变化规律及关键物候期RPI 特征分析

图7为第10号样本田时序RPI经过S-G滤波重构后拟合得到的RPI 时间序列曲线。从中可知，水稻的RPI 的数值变化范围为0.05～0.35。RPI 曲线变化趋势与水稻关键物候变化密切相关，有助于总结出水稻物候变化规律并识别出水稻关键物候期。

图6 第10号样本田时序RPI经S-G滤波处理前后对比图Fig.6 Comparison chart of the time-series RPI of sample field No.10 before and after S-G filter

1）出芽期至幼苗期。该时期的水稻秧苗矮小，分布较为稀疏，水稻下垫面为湿润的土壤。这也导致了体散射和二次散射分量较小，由平坦的土壤导致的面散射分量占优，但是整体的散射能量较低，RPI处于较低水平。

2）分蘖初期至分蘖中期。水稻分蘖期是决定穗数并为穗粒发育奠定基础的重要时期。水稻在这一生长阶段秧苗高度渐渐增加，稻叶的数量、叶长、叶宽和叶厚也均进行生长和分化，但此时的分蘖数以及水稻茎干长度还未达到峰值。由于水稻的分蘖和茎干的生长发育，植株自身结构之间的体散射逐渐增加，面散射分量逐渐降低，此时水稻冠层还未发育完全，冠层对散射能量的衰减作用较小，第一特征值λ1在此时近似达到最大值，使得RPI 在分蘖初期开始逐渐上升并在分蘖中期达到最大值。

3）拔节期至孕穗期。水稻分蘖末期与拔节初期在时间上有部分重叠。当水稻进入到拔节期，分蘖数目逐渐增加，稻叶发展完全，茎干长度逐渐伸长，形成更为茂密的冠层。尽管水稻冠层还未完全封闭，但对散射能量的衰减作用在逐渐变强。随着水稻生长至拔节中后期，水稻的分蘖数以及茎干长度达到最大值，水稻冠层有了进一步的发展，但还未达到最密。水稻生长至拔节后期至孕穗期，该时期属于营养生长和生殖生长并进阶段。孕穗时，幼穗包裹在叶鞘内，稻叶的形态并未弯曲。水稻茎秆已完全生长并达到峰值，此时水稻叶片在数量、长度和宽度方面均达到了最高点。从拔节期开始，由于水稻冠层的不断发育以及对微波能量衰减作用的不断增强，各个散射分量都迅速下降。在水稻孕穗期，尽管水稻冠层引起的衰减作用比下一时期（即抽穗扬花期）的要低，但由于水稻冠层已发育至一定规模，能量衰减的变化速率在孕穗期后逐渐变缓。反映到RPI曲线上，孕穗期是变化速率最快的点。

4）乳熟期至完熟期。乳熟期至完熟期为水稻生殖生长阶段。水稻生长至乳熟期之前，水稻还未开始灌浆，稻穗生物量小，稻穗形态直挺，水稻的株高达到最大值，形成了封闭的冠层。封闭的水稻冠层对能量的衰减作用达到最大，使得接收到的回波强度降到最低，λ1所代表的主要散射机制的强度也达到最低值，从而RPI 达到最小值。从乳熟期开始直至完熟期，水稻的生物量逐渐增加，稻穗形状弯曲，穗倾角逐渐增大，叶片开始干枯变黄，直至枯萎或脱落。水稻冠层的衰减作用相较于上一阶段减弱，各种散射机制的强度开始缓慢增加，3种散射机制中还是以体散射为主导，水稻冠层对于面散射和二次散射的衰减作用依旧很强。因此，在RPI 曲线中，RPI最小值所对应的水稻物候期为乳熟初期。

3.2.2 基于曲线导数法的水稻关键物候期识别

水稻关键物候期提取准则：首先根据RPI 最大值（RPI时序曲线一阶导数为0的点，且由正变负）识别出水稻分蘖中期所对应的DoY。其次根据RPI最小值（RPI时序曲线一阶导数为0的点，且由负变正）识别出水稻乳熟期DoY。然后确定RPI曲线的拐点（RPI时序曲线二阶导数为0的点，且由负变正），识别出水稻孕穗期的DoY。通常，粳稻的分蘖初期与分蘖中期的间隔时间为15～20 d，由此可得到水稻分蘖初期大致的DoY。分蘖初期之前即为水稻出芽期至幼苗期。RPI最大值与拐点之间为水稻拔节期至乳熟期，拐点与最小值之间为水稻抽穗扬花期，最小值之后为水稻乳熟至完熟期。

以10号水稻样本田为例进行关键物候期识别，最终识别出的水稻物候阶段见图7，包括出芽期至幼苗期、分蘖初期至分蘖中期、拔节期至孕穗期、抽穗扬花期、乳熟期至完熟期。

基于本文提出的上述水稻物候识别算法，对2015 年的19 个水稻样本田进行水稻关键物候期识别。

3.2.3 水稻关键物候期识别结果分析

3.2.3.1 评价指标

用xi表示识别出的水稻物候期所对应的DoY，用yi表示地面调查记录的水稻物候期所对应的DoY，以yi为真值，xi为待评价值，计算本文提出的用水稻物候识别算法识别出的物候期DoY 和地面同步实验物候期所对应DoY 的绝对误差（absolute error, AE）和均方根误差（root mean square error,RMSE），二者的表达式如下：

图7 粳稻生长的各个物候阶段（10号样本田）Fig.7 Various phenological stages of japonica rice growth(sample field No.10)

式中N为样本个数。

3.2.3.2 结果验证

图8 水稻样本田识别结果与野外观测物候比较Fig.8 Comparison of identification results of rice sample fields with field observation phenology

图9 19个水稻样本田获取的水稻物候反演结果与观测值对比分布图Fig.9 Distribution map of rice phenology retrieval results and observation values obtained from 19 rice sample fields

图8和图9是利用本文提出的算法识别出的水稻关键物候期DoY 与地面测量实验记录的DoY 比较结果。从中可以看出，识别出的分蘖中期、孕穗期和乳熟期的误差大部分在±8 d 以内，全部落在±16 d以内，在各个水稻样本田识别出的水稻关键物候期误差都处在一个较为合理的范围内，仅17号水稻样本田识别出的3个水稻关键物候期误差相对较大。由于水稻物候期所代表的是一个时间段，本文识别结果的误差处于合理的区间内，说明利用RPI曲线识别水稻物候期具有可行性。用RPI识别出的水稻关键物候期与观测值的均方根误差均小于6 d，其中乳熟期识别精度最高，为4.88 d，分蘖中期误差最大，为5.55 d，孕穗期为4.88 d。此外，研究区内粳稻的各关键物候期基本于同一时间到达，同一物候期的时间跨度在10 d 以内。杨振忠利用水稻全生育期的光谱特性结合K近邻、支持向量机等6 种机器学习方法，识别出分蘖期、拔节孕穗期、抽穗扬花期、灌浆成熟期4个生长发育阶段，整体识别精度在76%以上[31]；何泽基于多时相RADARSAT-2 数据，分别利用后向散射系数和Cloude-Pottier 等极化分解参数来建立水稻全生长周期反演模型，识别出了水稻移栽期、营养生长期、生殖生长期和成熟期，总体识别精度在82%以上[32]。与上述已有研究结果相比，总体来看，利用本文提出的RPI识别出的水稻关键物候期误差小，取得了令人满意的结果。

4 结语

本文利用多时相全极化SAR 数据对江苏省金湖县附近地区水稻种类进行识别，通过对水稻极化特征参数进行时序分析，创新性地构建出了能反映水稻物候变化特征的雷达物候指数（RPI），并利用其对水稻关键物候期进行反演。结果表明：利用香农熵可将籼稻和粳稻较好地识别出来，识别精度分别为92.38%和95.10%；利用本研究提出的雷达物候指数（RPI），反演出的水稻关键物候期日期与野外地面调查获得的日期相差全部在±16 d以内，得到了较准确的反演结果。本研究为深入剖析水稻雷达极化特征响应提供了可靠的理论依据。