运用Sentinel-1A时序数据和动态时间规整算法提取稻虾田空间分布信息

卢泽如 熊勤学 周雨顺 邓琪云 纪绍威 丁璐

摘要：针对稻虾田水稻移植期的不确定性给稻虾田空间分布信息提取带来的困难，运用SAR VH极化年时序数据、VV极化年时序数据以及调查区稻虾田各自极化的年时序数据通过标准化处理后，进行TWDTW方法（time-weighted dynamic time warping）计算，得到二波段（VH极化、VV极化）与调查区稻虾田年時序数据相似度空间分布栅格数据，再利用最大似然法进行监督分类，得到2016—2018年3年的监利县稻虾田空间分布栅格数据，分类结果与4块区域（龚场镇、福田寺镇、红城乡、白螺镇，面积分别是373.33、2 200.00、126.67、240.00 hm2）稻虾田实际分布图比较，其Kappa系数分别是0.79、0.80、0.75、0.90，介于0.75～0.90之间，说明运用Sentinel-1A时序数据和动态时间规整算法的提取方法是准确的。由于SAR传感器具有全天侯以及不受云层影响等特点，能定期获取研究区上空数据，相比光学传感器获得无云或者少云数据的不确定性，此方法能在实际操作中运用。

关键词：Sentinel-1A;动态时间规整算法;稻虾田空间分布信息;TWDTW算法;VH极化;VV极化

中图分类号：S127 文献标志码： A

文章编号：1002-1302（2020）18-0230-07

收稿日期：2019-10-25

基金项目：国家自然科学基金（编号：31871516）;公益性行业（农业）科研专项（编号：201203032）;湿地生态与农业利用教育部工程研究中心开放基金（编号：KF201701、KFT201906）;长江大学第十一批大学生创新创业训练计划（编号：2018288）。

作者简介：卢泽如（1999—），女，湖北黄冈人，主要研究方向为农业信息。E-mail：1284563860@qq.com。

通信作者：熊勤学，硕士，副教授，主要研究方向为农业遥感。E-mail：nxqx@tom.com。

近几年来，被农业农村部誉为“现代农业发展的成功典范”的稻虾共作生态种养高效模式在全国飞速发展[1]。根据全国水产技术推广总站发布的《2017中国小龙虾产业发展报告》，全国小龙虾（克氏原螯虾，Procambarus clarkii）养殖面积超过80万hm2[2]。如何研究稻虾共作模式的生态效应，为稻虾共作模式的可持续发展提供依据，是当前亟待解决的问题[3]。基于遥感影像数据提取稻虾田空间分布信息，被视为最便捷、最高效的研究手段，而目前对稻虾田空间分布信息的遥感提取研究还不多。魏妍冰等根据稻虾田独特的水域变化特征，利用自动水域提取指数（automated water extraction index）分别获取不同季节的水体，最后根据水体季相差异特征，用2017年8月和12月2景Landsat 8 OLI影像数据提取了2017年潜江市稻虾田的空间分布;这种冬季与夏季相比会增加的水域，视为稻虾共作田的方法有明显的异物同谱的错误，如低洼地带的中稻田和水生植被区[4]。造成提取稻虾田空间分布信息困难的原因主要有3个：一是同一地区稻虾田的水稻移植期差异很大。小龙虾市场价格因素和优良的水热气候条件使得水稻可生长的日期远大于实际生长期，水稻移植期为5月中旬至7月中旬，过长的水稻移植期导致不同稻虾田光谱特征差异很大，增加了遥感影像数据的提取难度;二是长重访周期和雨热同季的气候特征限制了光学卫星数据真实表达物候光谱特征的可能性[5];三是复杂的种植制度和粗光学遥感影像空间分辨率造成的混合像素严重制约信息提取的精度。

合成孔径雷达因不受光照和气候等条件限制，具备实现全天时、全天候对地观测的特点，特别是欧空局Sentinel-1A卫星C波段合成孔径雷达（SAR），实现了多极化（HV、HH极化）、高重访周期（12 d）、高空间分辨率[干涉宽幅（IM）模式：250 km，5 m×20 m分辨率]对地观测，海量数据加上免费政策，为稻虾田空间分布信息的提取提供了足够的基础数据。本研究意在以全国稻虾种养面积最大的县——湖北省监利县为研究对象，在充分分析稻虾田Sentinel-1A卫星C波段SAR数据时序特征的基础上，改进能够对时序数据进行动态时间归整（dynamic time warping，DTW）法，来消除长时间水稻移植期对不同稻虾田光谱特征差异的影响，高精度地提取监利县稻虾田空间分布信息。

1 研究区概况和数据预处理

1.1 研究区概况

监利县（图1）位于湖北省中南部，长江北岸，面积3 508 km2，四季分明，热量丰富，光照适宜，雨水充沛，雨热同季，无霜期长，是全国稻虾种养面积最大的县（2017年达到3.33万hm2）[2]，也是水稻种植面积最大的县（市）。

2017年8月，为分析稻虾田SAR数据的时序特征和验证结果的精度，在监利县的龚场镇、福田寺镇、红城乡、白螺镇分别选取4块验证区，并用GPS仪确定了每块稻虾田的边界（图1），4块区域稻虾田种植的面积分别是373.33、2 200.00、126.67、240.00 hm2。

1.2 数据的预处理

从欧洲航天局的网站（https：//scihub.copernicus.eu/dhus/#/home）上下载2016年1月至2018年12月涵盖监利县区域的Sentinel-1A卫星C波段SAR GRDH的格式数据，2景数据能涵盖监利县所有区域，共下载79 d共158景（2016年40景、2017年56景、2018年62景）数据，运用ESA Sentinel 1 Toolbox（Ver 1.1.1）软件做数据的预处理，数据预处理步骤如下：

（1）SAR数据合并：监利县跨越2景数据，因此在ESA Sentinel 1 Toolbox软件中打开同一天的2景GRDH数据，运用SAR MOSAIC功能，将2景数据合并;（2）SAR数据裁剪：加载SHAPEFILE格式的监利县行政区划电子数据，在Navigation窗口中缩放到合适位置后，利用spatial subset from view功能将涵盖监利县的数据裁剪下来;（3）后向散射系数计算：运用Radiometric下的Calibrate功能，选择output sigma band选项，计算VV、VH 2种极化、裁剪区域的后向散射系数;（4）图像降噪和滤波处理：为降低图像噪声水平，运用Multi-Looking功能进行多视处理，运用Speckle Filter功能，选取Refine Lee滤波方式进行滤波处理;（5）地形校正：运用Range Doppler Terrain Correction功能，选取SRTM 3Sec数字高程，坐标系统选取UTM WGS84 49N进行高程校正;（6）数据导出：将地形校正后的数据以ENVI格式导出，生成VV、VH 2波段数据;（7）多时序VV、VH极化数据生成：将2016—2018年数据按上述步骤进行处理，在ENVI软件中按年打开所有处理数据，运用Layer Stacking功能，将1年内所有VV极化或VH极化数据按日期顺序排列，生成一个多波段时序数据，波长用DOY表示（一年中的第n天），这样一年中有2个多波段时序遥感数据，一个是VV极化的，另一个是VH极化的。

1.3 分类结果验证方法

将结果数据和验证区调查空间数据运用GIS中的聚类分析生成只含0和1数字的栅格数据，0代表非稻虾区，1代表稻虾区，然后进行Kappa系数计算，公式[6]如下：

Kappa=Po-Pc1-Pc。

其中：Po=s/n。

Pc=（a1×b1+a0×b0）。

式中：n为栅格总像元;验证区栅格数据中为1的象元数为a1，为0的象元数为a0;结果模拟数据中为1的象元数为b1，为0的象元数为b0;2个栅格对应象元值相等的象元数为s。Kappa计算结果中通常Kappa值是在0～1之间，0.0～0.20表示极低一致性，0.21～0.40表示一般一致性，0.41～0.60表示中等一致性，0.61～0.80表示高度一致性，0.81～1.00表示几乎完全一致。

2 稻虾田提取原理与方法

2.1 稻虾田提取原理

影响地物后向散射系数的要素很多，有地物的粗糙度、土壤含水量、作物结构（生物量、叶面积系数、叶面、茎杆粗细等）、极化方式、频率和入射角等[7-8]。相同地物的后向散射系数变化规律基本一致，这是利用微波后向散射系数的时序特征进行地物分类的基本原理。水稻田后向散射系數的时序特征非常明显，即营养生长期（5—6月），因水稻淹水，后向散射系数偏低，与水体的时序特征相似，而进入生殖生长期（7—9月），其后向散射系数明显增加，到成熟早期达到最大[9];无论是稻虾共作还是连作方式种养的稻虾田，水稻收获后到下轮水稻移栽期（11月至次年4月），为了方便养殖小龙虾，稻田里一直会保持适当水位的水，因此这段时间的后向散射系数和水体一样，是偏低的。图2为2017年稻虾田、传统水稻田、水体、旱地、城市5种地物VH极化后向散射系数的年变化曲线。水体后向散射系数一年四季都很低（小于0.01），而城市后向散射系数一直偏高（大于0.12），旱田后向散射系数介于0.03～0.07之间，而传统稻田4—6月后向散射系数小于0.02，其他时间介于0.03～0.07之间。稻虾田后向散射系数的时序变化特征较其他地物有明显不同，其后向散射系数只会在6—9月偏高，其他时间都偏低。

2.2 稻虾田提取方法

尽管稻虾田后向散射系数的时序变化特征较其他地物有明显不同，但因为监利县光热资源丰富，大于20 ℃的气温持续133 d，属双季稻种植区，水稻可生长期大于实际生长期。而不同种植户因受龙虾养殖情况和市场因素影响，插秧期差异很大，从5月底至7月初都有可能插秧。2017年在监利县调查的4个地点，插秧期为5月23日至7月2日，每个点的HV极化后向散射系数的时序变化曲线有明显的平移变化（图3），而且最大值也有差异，因此常规的遥感分类方法不适用于稻虾田的空间分布信息提取。而DTW算法可以计算2个时间序列的相似度，尤其适用于不同长度、不同节奏的时间序列，能消除插秧期差异带来的后向散射系数时序变化曲线的平移影响。因此，本研究拟采用DTW方法提取稻虾田的空间分布信息。

2.2.1 DTW计算方法

DTW方法最初用来辨别2组语音的相似性，即计算2组时间序列的相似度，尤其适用于不同长度、不同节奏时间序列相似度的计算，近几年常用于遥感数据时序分析中。Petitjean等从理念的角度介绍了DTW方法，并说明了如何应用在遥感数据时序分析中[10]。Baumann等运用DTW方法消除了年际间物候差异对NDVI指数时序的影响，并成功对森林进行分类[11]。Gristina等利用DTW方法识别出Landsat NDVI时间序列中的异常值，并用时空插值方法代替，从而改进森林分类精度[12]。Maus等将时间限制条件引入DTW方法中，提出TWDTW方法（time-weighted DTW）来提高分类精度[13]。由此可知，DTW主要用来计算遥感数据归一化指数时序数据的相似度。为将DTW方法运用于SAR后向散射系数相似性，提出先将数据进行标准化处理，然后再进行TWDTW的方法，具体计算方法如下。

如计算2个点后向散射系数时序数据r=（r1，r2，…，rn）与 t=（t，t2，…，tm）之间的相似度。

（1）对后向散射系数时序数据进行标准化处理：

Ri=（ri-u）/σ。

u、σ分别为数组的r平均值和方差。

（2）用欧氏距离方法计算出数据R中的每个点与数据T中各点的距离，生成一个n行m列的距离矩阵D[n，m]。

（3）DTW相似度计算

DTW（i，j）=DTW（i-1，j-1）+min{D（i-1，j），D（i-1，j-1），D（i，j-1）}。

其中DTW（0，0）=0;

当i-j>K时，DTW（i，j）=0，K为时间限制，根据稻虾田中水稻移栽期不确定的特点，K值确定为50 d。

当i或者j等于n或者m时，DTW（i，j）为最终2组时间序列的相似度。

2.2.2 稻虾田提取计算流程 稻虾田提取计算具体流程见图4。

3 结果与分析

3.1 时序数据标准化处理对稻虾田空间分布信息提取的影响

图5为调查点2经过标准化处理与未标准化处理的TWDTW计算结果空间分布图，其值越低，表明2组数据的相似度越高，图上显示为黑色。从图5可以看出，未标准化处理的TWDTW计算结果（图5-B）与稻虾田空间分布差异很大，主要原因是影响SAR后向散射系数的因子很多，除了地物差异外，还有土壤水分、土壤类型、植被信息等，导致尽管其时序特征曲线非常相似，但不同区域的时序曲线中最高值、最低值差异很大，直接运用TWDTW计算，其他特征引起的差异高于时序特征引起的差异;而时序特征曲线经过标准化处理后，将时序特征曲线置于同一量级，能准确反映物候引起的时序差异（图5-A）。因此，对SAR时序数据作TWDTW处理之前，一定要进行标准化处理。

3.2 稻虾田空间分布信息提取结果验证

图6为4块区域（龚场镇、福田寺镇、红城乡、白螺镇，面积分别是373.33、2 200.00、126.67、240.00 hm2）稻虾田实际分布图和提取结果图。通过计算得出它们的Kappa系数分别是0.79、0.80、0.75、090，介于0.75～0.90之间，说明运用TWDTW计算后，再用最大似然法提取的结果还是能真实反映稻虾田空间分布的。从对比图中可以看出，其主要误差是边界模糊不清和稻虾田有中间镂空现象，这主要是由于SAR图像中存在着显著的相干斑噪声[14]，尽管图像进行过降噪和滤波处理，但它对提取结果的影响还是很大。

3.3 监利县稻虾田空间分布特征

监利县稻虾田种养面积飞速发展，2016年面积为1.88万hm2，2017年为3.22万hm2，至2018年为4.71万hm2，与监利县农业统计年鉴的统计数据（2016年种养面积为2.00万hm2、2017年为3.33万hm2、2018年为5.00万hm2）基本一致。图7为2016—2018年监利县稻虾田空间分布，可以看出，监利县稻虾田主要分布在洪湖西边的低洼地带（监利县的东边），沿长江民垸区，也就是由以前传统的水稻田和鱼池改建而来。2018年稻蝦田种养面积超过0.27万hm2的乡镇有黄歇口镇、白螺镇、棋盘乡、上车湾镇、尺八镇。2016—2018年各乡镇稻虾田种养面积统计见图8。

4 结论与讨论

本研究运用SAR VH极化年时序数据和VV极化年时序数据，与验证区稻虾田的相同极化年时序数据进行TWDTW计算后，得到二波段DTW数据，再用最大似然法成功提取2016—2018年监利县稻虾田空间分布信息，其Kappa系数介于 0.75～0.90之间，与实际状况基本一致。由于SAR传感器具有全天侯监测以及不受云层影响等特点，能定期获取研究区上空数据，相比光学传感器获得无云或者少云数据的不确定性，此方法可在实际操作中运用。由于数据来源于VV极化和VH极化数据，相比用单极化数据提取，其结果精度有一定的提高。TWDTW法能够对时间轴进行动态拉伸或压缩，对作物种植时间存在提前或延后的耕地时间序列也有较好的适应性，特别适合稻虾田这种水稻移栽期不确定的农田分类。本研究提取稻虾田空间分布信息的方法，克服了前言中提到的提取稻虾田空间分布信息的难点，不失为一种行之有效的方法。

不可否认，尽管SAR数据的空间分辨率和时间分辨率较以前有很大的提升，但运用SAR数据进行分类，与同级别光学数据的分类结果相比，精度偏低[15]。主要误差来源是SAR图像相干斑噪声，这会导致分类结果的边界模糊不清;而光学卫星数据尽管获取的概率低，但如果能获得无云数据，可以得到清晰的边界。因此，如果合理地运用两者之间的优点，则既能有清晰的边界，又能利用SAR定时获取数据的特点而得到高精度的分类结果，是提高稻虾田空间分布信息精度主要的研究方向。

随着目前光学卫星数据空间分辨率的提高，面向对象的分类法开始运用[16]。面向对象的分类法既能提取对象的光谱信息，也能对对象的纹理、几何特征、位置信息进行量化分类。因此，将光学卫星数据和SAR卫星数据进行融合，并结合TWDTW法进行面向对象的分类，是将来稻虾田分类信息研究的发展方向。

参考文献：

[1]秦尊文. 以“虾稻共作”模式为抓手推进体制机制创新——潜江市全国中小城市综合改革的观察与思考[J]. 中国发展，2016，160（6）：51-56.

[2]农业农村部渔业渔政管理局全国水产技术推广总站中国水产学会.中国小龙虾产业发展报告2019[J]. 中国水产，2019，7（9）：12-19.

[3]曹凑贵，江 洋，汪金平，等. 稻虾共作模式的“双刃性”及可持续发展策略[J]. 中国生态农业学报，2017，25（9）：1245-1253.

[4]魏妍冰，陆 苗，吴文斌. 基于水体季相差异的稻虾共作提取方法研究[J]. 中国农业资源与区划，2019，40（3）：14-20，34.

[5]Guan X D，Chong H，Liu G H，et al. Mapping rice cropping systems in Vietnam using an NDVI-based timeseries similarity measurement based on DTW distance[J]. Remote Sensing，2016，8（1）：19.

[6]Bradley A P. The use of the area under the ROC curve in the evaluation of machine learning algorithms[J]. Pattern Recognition，1997，30（7）：1145-1159.

[7]Graham A J，Harris R. Extracting biophysical parameters from remotely sensed radar data：a review of the water cloud model[J]. Progrom Physics Geography，2003，27（2）：217-229.

[8]Inoue Y，Kurosu T，Maeno H，et al. Season-long daily measurements of multifrequency （Ka，Ku，X，C，and L） and full-polarization backscatter signatures over paddy rice field and their relationship with biological variables[J]. Remote Sensing of Environment，2002，81（2/3）：194-204.

[9]Toan T L，Laur H，Mougin E，et al. Multitemporal and dual-polarization observations of agricultural vegetation covers by X-band SAR images[J]. IEEE Transactions on Geoscience and Remote Sensing，1989，27（6）：709-718.

[10]Petitjean F，Inglada J，Gancarski P. Satellite image time series analysis under time warping[J]. IEEE Transactions on Geoscience and Remote Sensing，2012，50（8）：3081-3095.

[11]Baumann M，Mutlu O，Richardson A D，et al. Phenology from landsat when data is scarce：using MODIS and dynamic time-warping to combine multi-year landsat imagery to derive annual phenology curves[J]. International Journal of Applied Earth Observation and Geoinformation，2017，54（2）：72-83.

[12]Cristina G，Joanne C W，Michael A W，et al. Integrated object-based spatiotemporal characterization of forest change from an annual time series of landsat image composites[J]. Canadian Journal of Remote Sensing，2015，41（4）：271-292.

[13]Maus V，Cmara G，Cartaxo R，et al. A Time-weighted dynamic time warping method for land-use and land-cover Mapping[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 2016，9（8）：3729-3739.

[14]Moreira A. Improved multi-look techniques applied to SAR and scan SAR images[J]. IEEE Trans on Geoscience and Remote Sensing，1991，23（4）：529-534.

[15]McNairn H，Shang J，Jiao X，et al. The contribution of ALOS PALSAR multi-polarization and polarimetric data to crop classification[J]. IEEE Geoscience and Remote Sensing，2009，47（12）：3981-3992.

[16]Ovidiu C，Mariana B，Gregory P A，et al. Object-based time-constrained dynamic time warping classification of crops using sentinel-2[J]. Remote Sensing，2019，11（10）：1257.