基于Sentinel-1A双极化时序数据的甘蔗株高反演方法

孙 盛 刘立露 胡忠文 余 旭

(1.广东工业大学计算机学院， 广州 510006; 2.深圳大学自然资源部大湾区地理环境监测重点实验室， 深圳 518000; 3.广东工业大学土木与交通工程学院， 广州 510006)

0 引言

甘蔗是我国主要的糖料作物，被农业部列为九大主要农作物之一，广东省作为甘蔗的主要生产区，2019年种植面积达1 030 065.01 hm2，占据全国的12.28%[1]。因此，多维度地监测甘蔗的生长变化规律，对广东甘蔗产业的可持续发展具有较好的推动作用[2]。甘蔗株高作为甘蔗生长情况、产量的重要评估标准之一，选择株高进行反演对于甘蔗的生长监测具有较好的代表意义。

传统反演作物参数主要通过实地调研，需要耗费较多人力物力与时间，遥感技术以其在时空尺度上的优势，逐渐在几种大宗粮食作物如水稻[3-4]、玉米[5-6]、小麦[7-8]、油菜[9-10]等得到了广泛应用，但应用于甘蔗的研究仍然较少。

目前，学术界较多采用光学遥感技术完成甘蔗种植信息的识别[11-12]、生物物理参数反演与估产[13]。但是，光学遥感图像受天气和气候影响较大，我国甘蔗种植多在中国南部广东、广西地区，这些地区空气湿度大，常年多云多雨，应用光学遥感技术进行甘蔗生长监测的效果较差。合成孔径雷达(Synthetic aperture radar, SAR)遥感[14]具有全天时、全天候的特点，可以较好地适用于广东地区农作物生长监测。近年来，基于SAR后向散射系数对甘蔗进行作物类型鉴定[15-16]和长势参数株高[17]等生长监测的研究有所增加。但是，后向散射系数过早饱和的特点限制了对甘蔗生长后期的监测，因此，需进一步挖掘SAR的极化特征在甘蔗株高反演中的潜力。

目前，已有不少利用极化特征开展其他作物的研究，如CHANG等[18]通过对全极化PALSAR(L波段)卫星数据进行极化特征分解，推得极化雷达植被指数(Polarimetric radar vegetation index,PRVI)定量反演灌木丛生物量，该研究认为灌木丛生物量与PRVI的相关性(R2=0.75)优于雷达植被指数(Radar vegetation index,RVI)[19]的相关性(R2=0.5)。WANG等[20]基于全极化无人机合成孔径雷达(Uninhabited aerial vehicle synthetic aperture radar，UAVSAR)数据的极化分解方法，为不同的散射机制定义SAR监测作物生长的植被生长指数。结果表明，散射特性可以随作物类型和物候发展阶段而变化，该植被生长指数与作物株高和生物量的地面测量值具有很好的相关性。尽管这些雷达植被指数可以很好地表征植被情况，但仅限于全极化或简缩极化SAR卫星数据；另一方面，这些研究都是基于经济成本较高的商业SAR卫星，不太适用于大范围区域作物监测的应用。

本文采用欧洲航天局(European Space Agency, ESA)提供的免费科研模式的双极化SAR Sentinel-1A时序数据，分析双极化雷达植被指数(DPRVI)[21]对于株高反演的可行性；分析讨论DPRVI对甘蔗不同生长期的响应，采用4种经典的经验回归模型定量反演甘蔗株高并对最佳反演模型进行验证。将DPRVI与3种经典的极化参数的时间模型进行性能评价，以验证本文提出的株高反演方法的效果和性能。

1 研究区域和实验数据

1.1 研究区域与研究作物概况

研究区域位于广东省广州市南沙区，中心坐标东经113°56′，北纬22°78′，年降雨量为1 623.6～1 899.8 mm，2017年、2018年与2019年甘蔗种植面积分别达到6 600.00、5 973.34、4 913.34 hm2，是该区域种植面积最大的农作物[22]。在广州市南沙区大岗镇庙贝村进行研究区域的数据采集，如图1所示。甘蔗生长周期一般分为发芽期、幼苗期、分蘖期、伸长期和成熟期5个时期。在2、3月播种，至当年11月或次年1、2月收获结束。其中分蘖期5月至6月上中旬，伸长期6月下旬至10月，成熟期11月至次年1、2月。因甘蔗伸长期较长，将6月下旬至8月下旬分为伸长前期，8月下旬至10月分为伸长后期。

图1 研究区域位置Fig.1 Location of study area

1.2 研究区域数据获取

通过多次调研获取了研究区域甘蔗生长周期内23景时间序列Sentinel-1A数据及卫星过境时间基本同步的3期实地数据采集；此外，在当地气候与农业气象中心的支持下，补充了3个调研时间节点之间的甘蔗生长数据。数据获取时间为2020年2月底(甘蔗播种期)至2020年12月底(甘蔗成熟期)。

1.2.1SAR数据

为避免不同传感器观测配置对观测结果的影响，选取Sentinel-1A卫星的IW模式一级产品(Level-1)中的单复视(Single looking complex, SLC)图像，其重返周期为12 d、空间分辨率为5 m(距离向)×20 m(方位向)，运行轨道为升轨右视；单次成像时，完整条带长约250 km，距离向上由3个子条带构成，每个子条带在方位向上有9个bursts并按方位向顺序排列(图2)。

图2 Sentinel-1A IW 图像Fig.2 Sentinel-1A IW image

1.2.2实地调研数据

Sentinel-1A时序数据覆盖了研究区甘蔗的整个周期，地面实地调研时间分别为与卫星过境采集基本同步的3个时间节点：2019年11月24日、2020年9月29日与2020年12月22日(表1)。为保证足够的纯净像元，减少其他作物的影响，即降低混杂像元，实验区域要求种植面积不小于20 hm2，并通过北斗高精度分米级移动平台Qmini A5(亚米级差分信号)采集经纬度等数据。每次调研获取实验区域12块代表性地块的甘蔗株高、伸长期、种植密度与现场天气情况等；每个代表性地块大小为2 m×2 m，以米尺测量地面至甘蔗自然状态下最高点的距离作为甘蔗株高，将测量得到的每个地块的株高平均值作为该地块的甘蔗株高。

表1 实地调研时间点及现场图Tab.1 Field research date and scene pictures

2 株高反演总体方案及双极化雷达植被指数

2.1 总体反演方案

首先对时间序列Sentinel-1A数据进行预处理，将SAR散射矩阵转换为协方差矩阵；再通过Cloude-Pottier目标分解方法[23]求得特征值与极化度(Degree of polarization, DoP)[24]，算出雷达双极化植被指数DPRVI；对甘蔗株高h进行反演，得出甘蔗生长最佳反演模型。对于模型采用决定系数R2与均方根误差RMSE进行性能评测，并将DPRVI与3种经典极化参数进行对比。总体设计方案见图3。

图3 株高反演总体方案Fig.3 Scheme of plant height inversion

2.2 时间序列Sentinel-1A 数据预处理

由于SAR卫星的成像体制，Sentinel-1A图像数据含有相干斑噪声[25]。因此首先需对Sentinel-1A原始图像进行预处理。

利用Python搭建欧洲航天局哨兵系列卫星产品的综合应用平台SNAP(Sentinel application platform)二次开发环境；结合SNAP平台中的子模块完成图像的预处理。

将图像从散射矩阵转换为2×2的协方差矩阵(C2)，该矩阵是一个二阶散射信息量，由双极化散射矩阵目标向量k[26]的空间平均生成，计算式为

k=[〈SVV〉 〈SVH〉]T

(1)

(2)

式中SVH、SVV——双极化数据VH、VV极化通道的散射矩阵

〈·〉——各向同质情况下随机散射介质的空间平均

上标*表示复共轭。

在距离向和方位向上以4∶1进行多视处理[27]，减弱相干斑噪声影响和减少后续数据处理量；采用refined Lee filter 方法进行极化滤波处理；最终，将这些像素地理编码到UTM投影坐标系中。

2.3 双极化雷达植被指数DPRVI的计算

2.3.1Cloude-Pottier目标分解与特征值归一化

(3)

其中

Λ=〈|SVV|2〉

α=〈|SVV|2〉/〈|SVH|2〉

由式(3)可知，协方差矩阵C2为半正定的Hermite矩阵，故Cloude-Pottier分解可表达为

(4)

式中ci——秩为1的独立矩阵

λi、ei——C2的实特征值与特征向量

从式(3)、(4)得到协方差矩阵的两个特征值为λ1=Λ,λ2=Λα。

由定义可知，特征值量化了主要的散射机制。λ1与VV极化通道相关，代表奇偶次散射强度；λ2与VH极化通道相关，代表多次散射(体散射)强度。当目标由奇偶次散射占主导时，有λ1≫λ2。将λ1归一化处理得到归一化参数P1，计算式为

(5)

2.3.2极化度(DoP)计算

BARAKAT[24]提出了DoP的表达式并将其用于全极化数据的NN协方差矩阵。针对双极化数据进行改进，计算式为

(6)

式中 Tr()——矩阵的迹操作符

根据极化度定义可知，它量化了不同散射机制之间的相对强度，当DoP=1时，认为目标处于完全极化状态，由单一散射机制占主导，当DoP=0时，认为目标处于完全去极化状态，由多次散射占主导，多种散射机制共存。

2.3.3DPRVI参数计算

随着作物冠层的发育，散射机制的类型也变得复杂。作物的发育初期，主要为土壤表面的漫散射(奇次散射类型之一)；发育晚期则由来自于冠层和土壤的多次散射占主导,散射类型的随机性增加。因此，从作物发育初期到晚期，DoP与P1均呈下降趋势。MANDAL等[21]发现这两个参数对作物生长状态的表征存在一定差异性，基于这两个参数得到一个新的参数，即双极化雷达植被指数(DPRVI)[8，28]，并成功用于反演油菜花与小麦的叶面积指数(Leaf area index,LAI)和植物含水率。另一方面，文献[28-32]均证明了叶面积指数与株高有较好的相关性，对于甘蔗，通过对比何亚娟等[13]的甘蔗不同生育期LAI分析与本文的实验结果，可发现甘蔗伸长前期与株高具有很好的正相关性，即株高越高，叶面积越大。但随着后期叶片衰老、叶面积减少，株高并不降低甚至在增高，二者的正相关性不明显。综上，本文将DPRVI用于甘蔗伸长前期状态的监测并反演株高；同时对该参数在甘蔗伸长后期的株高反演进行初步分析和讨论。DPRVI的计算式为

DPRVI=1-DoP·P1

(7)

从定义可知，DPRVI∈[0,1]，DoP和P1的乘积对应于奇偶次散射的比例，通过单位1减去该比例可以推得目标生长过程的多次散射变化，即散射随机程度。理论上，当目标完全随机散射时，DoP=0(完全去极化)，λ1=λ2，即P1=0.5，此时DPRVI=1，等价于作物冠层完全发育状态下的散射特征；在光滑的裸露表面，散射现象主要为布拉格散射时，有DoP=1(完全极化)，λ1=λ1+λ2，λ2=0。此时DPRVI=0，等价于裸土的散射特征；数值在0～1之间的变化反映了作物生长的过程。

3 反演方法实现及实验分析

3.1 双极化雷达植被指数(DPRVI)对甘蔗株高(h)的反演

通过实地数据采集和咨询当地气候与农业气象中心及蔗农，得到不同生长期甘蔗平均株高，见图4。分析发现，株高随生长时间变化明显，通过比较，以Sigmoid函数进行拟合效果最佳，决定系数为0.980。从拟合曲线可看出，从发芽开始到伸长期前期甘蔗呈加速生长趋势，株高迅速增加，这是因为此期间光合作用逐渐增强，在此阶段，甘蔗生长特点为叶面积迅速增大，节间伸长增粗；进入伸长后期，甘蔗拔节变缓，株高增加变缓，这是由于此期间甘蔗吸收的营养转向蔗茎的增粗与糖分的积累。

图4 甘蔗不同生长期株高的Sigmoid函数拟合模型Fig.4 Sigmoid function fitting model for plant height of sugarcane at different growth stages

基于2—12月的时间序列Sentinel-1A数据，将株高拟合曲线与对应研究区域3个甘蔗田(甘蔗田1、甘蔗田2、甘蔗田3)的植被指数DPRVI绘制于图5。其中，甘蔗田1、甘蔗田2与甘蔗田3的种植面积分别约为3.33、4.00、3.67 hm2。通过实地测量，3块甘蔗田的平均行距约为100 cm，平均株距约为30 cm。

图5 甘蔗h与DPRVI随播期的动态变化曲线Fig.5 Dynamic changes of h and DPRVI with sowing date of sugarcane

由图5可知，在2、3月DPRVI的值较低，这是由于南沙甘蔗播种采取覆膜技术，田地经整土覆膜后以有序的奇偶次散射λ1为主导；随后，甘蔗进入幼苗期，体散射分量λ2增加，DoP降低，但仍以奇偶次散射为主导(λ1≫λ2)，DPRVI从3月到4月下旬的缓慢上升也验证了这一变化；而后DPRVI加速上升，于6月中旬达到峰值，这时甘蔗到了分蘖期，分蘖加快，冠层的快速发育使得体散射分量λ2进一步增加，此时田地表现为多种散射类型并存，散射随机程度最高；之后DPRVI下降，是由于甘蔗进入伸长前期，甘蔗分蘖变慢，节间加速分节与增粗，叶面积迅速增大，使得奇偶次散射分量λ1增加；到8月下旬后，DPRVI缓慢上升则是由于甘蔗进入伸长后期，分节与增粗变缓，叶片数积累增多，甘蔗冠层逐渐完全发育，体散射分量λ2增加，奇偶次散射分量λ1减少，DoP变小；进入成熟期后，DPRVI或因甘蔗叶枯萎掉落减低体散射分量λ2而有所降低，数值不再有较大波动。综上，图5中DPRVI的变化规律较好地反演了甘蔗的物候周期。在甘蔗生长期的分蘖期之前与伸长后期之后两者呈正相关关系，伸长前期两者呈负相关关系。鉴于以上分析，采用分段函数与4种经典经验回归模型(线性、二次多项式、指数、对数)针对不同生育期进行株高回归模型反演，以决定系数R2和均方根误差(RMSE)评价反演性能，结果如表2所示。

由表2可以看出，DPRVI与株高有较显著相关性，综合考量计算量、R2与RMSE性能，以二次函数模型反演效果最好，决定系数分别为0.882、0.625与0.716；将不同生育期的回归曲线绘制如图6所示，可以看出，幼苗期、分蘖期的相关性最佳，其次为伸长后期与伸长前期。

为进一步验证模型的可靠性，使用甘蔗田4(面积大约为20 hm2，平均行距约为100 cm，平均株距约为30 cm，位置见图1)作为验证数据，选择预测效果最佳的幼苗期、分蘖期二次函数模型进行比较(图7)。可以看出，相对于甘蔗田1、甘蔗田2和甘蔗田3，甘蔗田4的数据离散程度略大，这是由于甘蔗田4较少的耕种面积，受到了更多混杂像元的影响。接着进行预测值与实测值的比较(图8)，求得决定系数R2为0.839，平均绝对偏差为7.4%，说明利用该模型依然可以较好地反演甘蔗株高与DPRVI的关系，基本满足DPRVI对伸长前期(分蘖期之前)的甘蔗株高监测的动态需求，保证利用双极化雷达植被指数预测甘蔗伸长前期株高的准确性与可靠性。

表2 甘蔗各生长期h与DPRVI的关系Tab.2 Relationship between h and DPRVI in sugarcane growth period

3.2 不同极化参数对比分析

绘制2—12月的时间序列Sentinel-1A数据，求得3个甘蔗田研究区域的DPRVI、交叉极化与同极化通道的后向散射系数强度(σVH°、σVV°)、不同极化比值(σVH°/σVV°)、雷达植被指数(RVI)的时间模型，如图9所示。

由图9可知，3种经典的极化参数σVH°、σVV°、σVH°/σVV°与RVI的时间模型在3个甘蔗田的伸长前期均表现出一定的生长响应。其中σVH°、σVV°与h呈正相关性，σVH°/σVV°、RVI与h呈负相关性。但是，从甘蔗进入伸长期开始，该3种参数均已经趋向稳定，难以反演伸长期之后株高的生长变化。这是由于甘蔗生长早期散射机制相对简单，而到了后期，由于随着甘蔗冠层的发育分蘖，散射机制复杂性增加，后向散射系数达到饱和，基于后向散射系数的σVH°、σVV°、σVH°/σVV°与RVI也无法较好地表征株高，而基于极化特征求得的DPRVI在甘蔗不同生育期仍然保持了与株高较好的相关性。因此，在反演甘蔗作物时，DPRVI比Sentinel-1A后向散射系数及其推导得出的参数性能更加优异，对于指导甘蔗生产具有一定应用价值。

图6 甘蔗不同生育期株高最佳反演模型Fig.6 Optimal inversion model of plant height in different growth stages of sugarcane

图7 伸长前期株高反演模型与实测甘蔗株高的关系Fig.7 Relationship between inversed plant height at early growth stage and measured plant height of sugarcane

图8 伸长前期甘蔗株高实测值与预测值比较Fig.8 Comparison of measured and predicted height of sugarcane plant in early growth period

图9 不同甘蔗田的极化参数(DPRVI、σVH°、σVV°、σVH°/σVV°、RVI)时间模型Fig.9 Comparison of polarization parameter (DPRVI,σVH°,σVV°,σVH°/σVV° and RVI) time models in different sugarcane fields

4 结论

(1)将时序数据Sentinel-1A求得的DPRVI与3种经典极化参数进行对比。结果显示DPRVI能较好地反演甘蔗生长前期的株高变化。

(2)采用4种经典的经验回归模型对甘蔗株高反演建模。结果显示分蘖期之前两者相关性最高，二次多项式模型拟合效果最好，决定系数R2达0.882，均方根误差达0.118 cm。

(3)使用当地气候与农业气象中心数据进行对比分析，决定系数R2达0.839，平均绝对偏差为7.4%，验证了Sentinel-1A时序数据对粤港澳大湾区的甘蔗作物的连续监测是有效的，对于大规模、低成本的甘蔗作物生长状态监测，Sentinel-1A数据是一种较好的数据源。