地形约束下的地块级水稻分层提取方法

李梦秋杨树文骆剑承石含宁付昱凯

(1.兰州交通大学测绘与地理信息学院，甘肃 兰州 730070；2.地理国情监测技术应用国家地方联合工程研究中心，甘肃 兰州 730070；3.中国科学院空天信息创新研究院/遥感科学国家重点实验室，北京 100101；4.中国科学院大学，北京 100049；5.中国地质大学(武汉)计算机学院，湖北 武汉 430074)

引言

随着光学遥感技术在水稻识别中的应用日渐成熟[1]，多时相遥感数据能够进一步提升水稻识别精度[2]。利用光学光谱反射率、植被指数和水稻的物候期可以实现水稻制图，Xiao等提出，在水稻生育内期可以利用NDVI(Normalized Difference Vegetation Index，归一化植被指数)和NDWI(Normalized Difference Water Index，归一化水指数)来识别水稻田；随后其又利用LSWI(Land Surface Water Index，地表水提指数)加一些阈值等于或大于NDVI或其他增强型植被指数(LSWI-EVI或LSWI-NDVI)突出水稻生育期内水信号，设置适宜本地水稻生长条件的阈值来识别水稻[3，4]。上述方法表明了植被指数在水稻不同生长期内的变化能够明显将水稻区别于其他作物[5，6]，但受天气影响，在多云多雨的西南山地区，水稻识别在最佳物候期无法保证有高质量成像，使得该地区光学遥感时间序列难以构建，仅用光学遥感不能满足高精度水稻识别[2，7-9]。

与光学遥感相比，微波遥感不受光照、云雾等天气条件的影响，全天时、全天候检测，并且有穿透云层反映地面目标的电磁和结构特性的能力。多时相SAR数据受天气影响较小，可以完整捕捉水稻生长期内的物候信息，进一步提高水稻识别的精确性[10，11]。SAR数据在极端天气地区的水稻识别中表现出极大潜力[12-15]，但大多水稻识别研究主要针对具有大面积且连续的水稻种植模式的平原地区，而在西南地区水稻种植大多在山地、丘陵地区，水稻种植地块细小破碎且分布零散[16，17]，空间连续性不足，大尺度对象无法满足内部对象一致性，常伴有混合像元，导致水稻识别结果经常出现错分和漏分现象。Xu等构建了一种通过NDVI、NDWI以及VH极化运算的水稻识别模型[18]，有效地将光学和SAR影像结合使用，但未能重点关注到西南地区因地形导致水稻种植离散度较高，物候期在不同海拔具有差异性[19]。因此，在西南山地区域开展地块级水稻制图方法研究具有重要的意义和价值。

综上所述，本文利用水稻移栽时因蓄水造成地表粗糙度降低，SAR后向散射系数呈现下降再上升的“V”字趋势将Sentinel-1VH时序的后向散射系数结合Sentinel-2的NDVI、NDWI时间序列计算出每个耕地地块的水稻指数(SPRI)，在此基础上加入高程数据作为约束条件，研究区内耕地地块按照水稻生育周期对应的高程划分，分别划分为0～400m水稻地块及400m以上水稻地块，根据样本数据在不同高程分层下进行多阈值二进制分类。该方法针对我国西南地区地形复杂可以有效减少研究区内因地形引起的水稻识别错误从而提高了分类精度。

1 研究区概况

1.1 研究区位置

双河镇隶属于重庆市荣昌区南部，东西最大距离10.0km，南北最大距离9.5km，总面积88.6km2，见图1。受亚热带季风气候影响，全年平均气温17.9℃，年均降水1111.6mm，年平均日照时数1282h。双河镇地处川中南丘陵区，地形地貌以浅丘为主，地势大致为北低南高，东高西低；地块破碎分布零星，农作物种植结构复杂多样。研究区受四川盆地和云贵高原气候的相互影响，光热同季，气候温和，雨量充沛，是水稻主要生长区。

图1 研究区地理位置

1.2 数据源及预处理

1.2.1 Sentinel-1/2数据及预处理

本文选用的Sentinel-1数据源是由阿拉斯加卫星数据中心(https：//search.asf.alaska.edu/)提供的一级产品Sentinel-1 A GRDH地距多视影像，实验选取VH极化方式，入射角为38.278°。影像时间覆盖范围为2022年1—9月、11月，共17期雷达数据。为了进一步配合SAR数据进行水稻制图，本次实验还通过GEE(Google Earth Engine)平台获取了双河镇哨兵二号Level-2A数据，其分辨率为10m，并在该平台上通过对蓝、绿、红和近红外波段计算得到了NDVI和NDWI 2个光谱指数。

Sentinel-1影像在SNAP软件中预处理，主要包括6个步骤：轨道矫正，热噪声去除，辐射定标，地形矫正，噪声去除，以及分贝化。

1.2.2 其他辅助数据

1.2.2.1 其他遥感数据

基于重庆市农业产业数字化底图项目支持，由其项目组提供GF-2卫星遥感影像，其空间分辨率为0.8m。在GF-2卫星遥感影像中可清晰获取研究区道路、建筑、耕地、水域等地物边界，为后续耕地地块提取及样本补充提供了精细的空间特征。

1.2.2.2 地块数据

基于GF-2遥感影像的地块提取，将耕地形态进行分区，分为边缘和纹理特征；分别采用HED(Holistically-Nested Edge Detec-tion)模型和D-LinkNet模型对规则耕地区和坡耕地区进行地块提取[20]，得到研究区耕地地块分布信息。

1.2.2.3 地面调查数据

当地水稻作物物候信息由重庆市农委局提供作为参考，见表1，样本数据采用2022年3月和6月开展的作物类型实地调研结果及根据GF-2、Google Earth影像配合土地利用数据进行的样本扩充，共计获取水稻样本1123个。

表1 重庆市渝西地区水稻物候信息

1.2.2.4 高程数据

为引入高程数据作为研究区水稻分类的约束条件，本文选择的DEM数据为修正后的ASTER GDEM V2，其空间分辨率为30m，全域高程起伏范围介于266～597m，相对海拔较高区域多集中于南部。

2 研究方法

在探究地形对水稻物候期的影响和以耕地地块作为最小识别单元来加强水稻指数对水稻生育期内各特征的量化的研究分析基础上，本文构建了基于地形约束的地块级水稻分层提取方法，其具体流程如图2所示。

图2 基于地形约束的水稻分层提取流程

本次实验的研究方法主要包括4大部分。地块级光谱时序构建，将处理所得到的VH、NDVI、NDWI值赋到对应地块上，得到地块尺度的光谱特征；水稻指数计算，计算得到地块级水稻指数(SPRI)，SPRI∈(0，1)，值越趋近于1，是水稻的可能性越大；高程分层约束，修正研究区矢量边界保证边界地块完整性，配合当地水稻物候信息，带入移栽-高程公式得到最佳高程分界值，以高程作为约束条件将研究区划分再分别提取水稻；多阈值水稻分类及精度验证，在不同高程内设置不同分类阈值，在得到的结果中选取样方进行精度验证。

2.1 地块级光谱时序构建

在实现复杂地表条件下的精准农业过程中，面临的一大挑战是大多数遥感影像分类与提取主要基于像素单元，这种方法往往忽略了同一地块内临近像元之间的相关性。这导致了形态位置与实际耕地地块的不完全对应，从而削弱了地块内同种作物的同质性和地块间不同作物的异质性。针对这一问题，地块级光谱时序构建的方法显得尤为重要。根据Sentinel-2/Sentinel-1多时相数据，本文以地块为单位构建NDVI、NDWI光学特征及VH后向散射系数曲线。

以地块作为最小单元不仅可以整合光学影像所获取的水稻空间特征，还能与SAR影像所捕捉的水稻地块上的物候特征相结合[3]，更清晰直观地反映水稻的位置和空间形态。SAR的VH极化方式曲线随着水稻的生长周期呈现“V”字形变化，在水稻播种期和生长初期由于地块以蓄水为主，导致后向散射系数的反射源来自水面而非植被，此时数值较低。随着水稻进入生长期，植被体积和密度增加，此时后向散射系数逐渐增加。水稻进入抽穗期至成熟期时植株达到最大体积和最高水分含量，此时数值达到峰值，但随着植株逐渐干燥，系数会略有下降。在收获期，由于水稻被收割，制备大量减少，后向散射系数显著下降。水稻的NDVI和NDWI有效反映了植被的覆盖程度和水分含量，其变化趋势如图3所示，在播种和初生期，由于植被覆盖度低，NDVI和NDWI值都相对较低。进入生长期，随着植被的发展，这2个指数值逐渐升高。在抽穗至成熟期，NDVI值达到峰值，而NDWI值在这段时间内保持较高或略有下降。在收获期，随着植被被收割，NDVI和NDWI值都显著下降。同时，通过建立地块级的时序特征，能够加强同一高程下水稻地块的同质性，同时突出不同高程下水稻地块间的异质性。

图3 水稻地块后向散射系数曲线

图4 水稻地块光谱指数曲线

2.2 地块水稻指数(SPRI)构建

农作物的光谱曲线具有明显的周期性变化规律，后向散射强度对应的水-土-植被组成的时间变化是水稻识别的关键特征[18]。水稻在生长周期内由土壤、水、植被组成的动态混合物，而SAR后向散射机制主要受土壤水分、植株盖度的影响。在淹水期，由于水稻面积小，覆盖稀疏，此时后向散射强度主要受地表水的影响，因此水稻移栽期的后向散射值远低于未灌溉的其他植被或作物，水稻成熟期时后向散射强度与其他植被无显著差异。本文以地块为基础构建水稻指数，完整地捕捉不同高程内水稻地块上各指数完整的变化趋势，有效消弱“同谱异物”与“同物异谱”带来的分类错误。

SPRI是基于Sentinel-1 VH后向散射对不同土地覆盖类型时间变化特征的分析，根据水稻的生长周期特征，图5反映了在水稻的淹水阶段后向散射值低且接近水的后向散射值；水稻生长阶段的后向散射值逐渐升高且大于其他作物最高值接近植被的后向散射值。该趋势表明，在水稻的淹水阶段，其后向散射系数受蓄水影响较大，数值远远小于其他植被或作物，在成熟期时其后向散射系数与别的植被或作物又无明显差异。因此，利用水稻后向散射系数在水-土壤-植被组成的时间变化可以有效区别其它作物。SPRI对这些特征进行了量化[18]，其计算方式：

图5 水稻在生长周期内Sentinel-1后向散射系数时序变化特征

SPRI=f(D)×f(W)×f(V)

(1)

式中，f(D)为扩大水稻和其他作物之间的差异，用sigmoid函数计算了后向散射范围D=(P2-P1)和V线与W线的深度(v-w)之间的关系；f(W)用于计算p1和W线的接近程度；f(V)用于计算p2和V线的接近程度。f(D)、f(V)、f(W)的取值范围均为0～1。

(2)

(3)

(4)

式中，v和w分别表示V线和W线的后向散射值；p1和p2分别是生长期间的局部最小和最大强度。D、v、w、p1、p2的单位为dB。

2.3 行政区划与种植分区

我国水稻种植区非常广阔，各地自然生态条件复杂多样，社会经济条件各不相同，因此，水稻在不同区域的生长条件差异较大。行政区划在大范围上遵循了地域分异规律，按照人类活动特点和区域经济特点一定程度保证了不同行政区内各自的特色种植制度，可以更好地因地制宜进行农业生产布局。种植分区是在行政区划的基础上进一步参照地形、地势、地貌等本底信息来划分区域，缩小了各种植区内部作物类型、种植物候差异。

对照任务区内的省级、市级、区县级和乡镇级的矢量边界和高分2号影像来对乡镇边界进行微调，保证在边界耕地地块的完整性。

水稻的种植分布同时受到自然条件和人为管理的共同约束，早在1983年龚高法等[22]就提出了一个地区作物发育期的早晚受到当地气候、植被、土壤、地形等条件的影响。由此可知，一个地方的大气候条件制约于纬度、经度和海拔高度的共同作用。在此基础上王琛智等[19]研究得出，地形因素直接影响着水稻生长环境的温度、光照时长的空间分布，其中在地形因素中高程对水稻生长的影响最为突出，并给出水稻移栽期计算公式：

y移栽=-0.023×x高程+36.835

(5)

式中，y移栽为水稻移栽日期的DOY数值；x高程为该区域的高程值，m。根据式(5)和重庆市农技站提供的渝西地区水稻物候信息，见表1，得知双河镇不同种植区水稻移栽期前后时间差为10d，可计算得出400m的高程差是物候期的一个差异带。因此，本文在乡镇的行政分区的基础上进一步根据地形、地貌等自然特征进行二次分区，依据高程信息将研究区内划分为高程在0～400m的种植区及400m以上种植区，保证各种植区内地理条件的趋同性。

2.4 多阈值分类

提取出的耕地地块上均赋有土地利用类型、高程、水稻指数等属性，根据研究区当地水稻种植物候条件得知，水稻在海拔为0～400m时，SOS1(水稻最早的返青或出苗期)为95DOY，SOS2(水稻最晚的返青或出苗期)为110DOY，EOS1(水稻最早衰老或收获期)为263DOY，EOS2(水稻最晚衰老或收获期)为268DOY；水稻在海拔为400m以上时，SOS1(水稻最早的返青或出苗期)为105DOY，SOS2(水稻最晚的返青或出苗期)为120DOY，EOS1(水稻最早衰老或收获期)为263DOY，EOS2(水稻最晚衰老或收获期)为268DOY。依据高程差产生的同种水稻物候差异结合式(6)，将研究区的地块按照高程划分为0～400m及400m以上2个子区域。

根据土地利用迁移数据，将类型为耕地的地块筛选出来与实地采样点数据和目视解译结果取交，由于SPRI值表示地块种植水稻的概率，因此，可以用一个适当的阈值将每个地块划分为是水稻地块和非水稻地块。对不同高程区域内的SPRI设置阈值，根据水稻样本数数据得到水稻在0～400m海拔的地区SPRI≥0.5031，在400m以上海拔的地区SPRI≥0.5483；依照对不同海拔的SPRI值进行多阈值水稻分类，将各部分水稻分类结果合并作为整个研究区的水稻分类结果。

3 结果与分析

3.1 耕地地块提取结果及分析

利用高分2号遥感影像基于深度学习提取的双河镇耕地地块共49977个，见图6。双河镇的耕地地块在空间分布上呈现平缓地区地块集中且密集，海拔起伏较大地区地块稀少。海拔400m及以上有8492个地块，仅占总耕地地块面积的17%；其余地块均分布在0～400m海拔内。海拔较低区域的耕地大多连续分布，地块面积较大且形态较相似，而高海拔地区的地块分布零星且细小破碎，二者在面积、形态上差别较大。从地块提取的细节可以看出，见图6，地块边界明显，与实际地物符合度较高，内部结构相对统一，总体提取效果良好，以此为基准提取作物能够有效改善边界混合像元引起的错分、漏分现象。

图6 耕地地块提取结果

3.2 水稻分层提取结果及分析

本文提取的重庆市荣昌区双河镇2022年水稻种植情况结果如图7所示。双河镇水稻种植分布比较广泛，但主要集中在海拔0～400m较平缓地区，水稻地块共有34366个；随着海拔的升高，水稻种植地块数量逐渐减少，在海拔400m以上地区存在较少水稻地块且分布比较零星，水稻地块共有7401个。

图7 双河镇水稻地块空间分布

双河镇地势高低起伏，山地和平原相互交错，地形复杂多样。全域地势大致为南高北低，东高西低。从双河镇提取结果的水稻空间分布来看，海拔0～400m的区域内，水稻占耕地的比重约为38.9%，是水稻的主要种植区域，该区域地形平坦、水资源丰富，地块蓄水良好并能保持充足的湿度，是水稻灌溉、排水和施肥的理想地形；在海拔400m以上的区域，由于地形起伏较大，地块蓄水困难，容易流失水分和肥力，加之海拔较高土地开垦困难，鲜有农户在此开展作业，因此海拔400m以上区域水稻种植分布零星且稀碎，水稻占耕地的比重约为2.0%。

3.3 提取精度验证

本文将研究区按照最大范围边界分成了119个1000m×1000m的规则网格，并在研究范围内随机选取了4个网格作为样方展示水稻分类结果，见图8，其中海拔在0～400m的有2个样方，分别为1号样方，见图8a，2号样方，见图8b；海拔在400m以上的有2个样方，分别为3号样方，见图8c，4号样方，见图8d。

图9第2列展示了图8a～d 4个样方的SPRI值，高SPRI值(蓝色)与GF-2影像中的淹水区域(图9第1、3、5列的真彩色GF-2影像中的暗区)相吻合，这些被淹没的区域对应了水稻种植的淹水阶段；低SPRI值(红色)与其他植被覆盖区域相符。

图9 4个样方的高分影像和SPRI特征图

结合目视解译和实地样点的样方内11912个参考地块作为验证样本，采用准确率(Accuracy)、精确率(Precision)和召回率(Recall)和F1_score对样方内的水稻制图结果进行验证，见表2。从精度验证表可知，海拔在0～400m的F1_score达到92.1%，海拔在400m以上的F1_score达到90.8%，分类效果较好，见图9。

表2 水稻提取精度验证F1_score

本文以海拔分层为统计单元，由于2022年的统计数据存在滞后的情况，根据重庆农委给出的2021年水稻种植面积统计数据和2023年1月重庆市荣昌区发展和改革委员会发布的2022年粮食作物增幅比重数据，计算得到2022年水稻种植面积统计数据大约为218km2，通过双河镇与荣昌区全域面积比计算得到双河镇2022年水稻种植面积统计数据大约为17.7km2。根据本文方法分区域提取结果为海拔在0～400m地区水稻种植面积为20.1km2，海拔在400m以上地区水稻种植面积为3.58km2，合计全域水稻种植面积为23.68km2。由于本文是在耕地地块提取结果的基础上进行的水稻指数计算和分类，尤其是在地形起伏多变的地区有较多细碎且分布零星的地块不易被提取，导致水稻制图面积和统计面积存在差异同时也要结合当地种植模式的调整，荣昌区双河镇近年逐渐增加了“稻虾供作”模式，稻田蓄水时间会久于普通稻田，这也会让水稻提取结果与实际情况有出入。

4 结论与展望

为了验证本文的分层级水稻分类方法的有效性，对比无高程约束条件下进行单阈值分类。通过对比高程约束下多阈值水稻分类方法和无约束单阈值水稻分类方的结果精度和分类效果来评估不同分类方法的优缺点，采用样方验证和统计数据进行对比分析，基于实地采样数据和目视解译结果对2种方法在双河镇的制图结果进行验证。

由图10对比可以看出，在高程约束下的水稻地块提取细节上与验证样本更符合，无约束单阈值的错提主要因为在实际种植中不能避免在>400m区域存在一些自然植被与<400m区域水稻有相同物候期，从而增加了错提情况。因此，在不同高程下设置不同阈值进行水稻提取有效的针对物候期进行了划分，更适合因高程变化所导致的水稻物候期差异明显的西南地区，减少了水稻错提的现象。

图10 高程约束多阈值水稻地块提取与无约束单阈值水稻地块提取结果比较

无约束单阈值提取结果精度验证表明，见表3，高程约束下多阈值分类精度和无条件约束单阈值分类精度F1_score最高分别为92.1%和89.2%，且各个样方内精度也是约束下多阈值分类精度高于无约束单阈值分类精度。

表3 高程约束下多阈值分类与无条件约束单阈值分类精度验证对比

西南地区地形起伏多变导致不同海拔之间的水稻物候期产生差异，故无法在双河镇全域用同一个标准进行水稻分类。根据地理学第一定律空间相关性，即通常情况下两地物之间的距离越近那么相关性越大，本文按照高程分层使得每层内的自然条件的属性具有相对一致性，在相同海拔内水稻所受到的气候条件相同，在生长过程中生长周期一致；根据地理学第二定律空间异质性，即多个区域之间存在差异，在不同海拔间高海拔地区受小气候影响较为严重导致同一作物在生长周期上与平原地区不同。因此，按照高程将研究区分层是将复杂问题解构成简单问题的过程，保证了每一区域内的条件和影响因素具有相对一致性，减少了同物异谱的错分和漏分情况，使得研究区水稻分类精度提高。

高程约束下多阈值分类整体而言，海拔在0～400m的分类精度要高于海拔在400m以上的分类精度，F1_score分别为92.1%和90.8%。这是由于在耕地地块提取过程中平原地区耕地特点更规则且连片，而在山地地区种植条件有限，会出现分布零星且不规则的耕地地块进而导致0～400m区域耕地提取精度高于400m以上区域。加之海拔0～400m地区的气候相对一致和稳定，而400m以上地区受局部小气候影响，气候条件相对不太稳定致使雷达信号在接收时大气损耗较大，后向散射系数变化趋势较不稳定，因此，在低海拔地区的水稻提取精度比较高海拔区域的提取精度要高。

本研究针对西南地区存在的地形复杂、地块破碎的问题，引入高程作为研究区的约束条件对不同高程区域内地块进行多阈值水稻分类，通过对比无约束单阈值水稻分类结果，得到以下相应结论。

在地形复杂地区的分类中，“同物异谱”现象更加严重，引入高程约束条件，将复杂场景简化，减弱了同一高程区域内的差异性，从而减少了“同物异谱”导致的错分漏分问题，提高了分类精度。

不同高程区域的水稻生长环境不同导致SPRI值存在差异，多阈值分类进一步强调了具体问题具体分析，根据不同海拔下后向散射系数的变化设置适合该区域内的水稻阈值分类条件。

利用高程作为分层的约束条件虽然有效提高了研究区内的水稻提取精度，但在实际情况中还有许多因子可以作为约束条件去考虑，如坡度、坡向以及水热条件，在不同的实际影响情况下选择不同的因子作为约束条件或几个因子组合作为约束条件是本次研究未深入的地方，同时雷达的后向散射系数与这类因子之间的变化关系都是在后续学习中会持续深入探讨的，进一步提高分层约束在西南地区的可实用性。