基于Sentinel-1时间序列的河冰表面物候提取研究——以黄河内蒙古段为例

张 鑫，牟献友，翟涌光，张宝森，冀鸿兰

（1.内蒙古农业大学水利与土木建筑工程学院，内蒙古 呼和浩特 010018；2.黄河水利委员会黄河水利科学研究院，河南 郑州 450003）

0 引言

冰凌洪水灾害是高纬度河流特有的一种易发自然灾害，河冰物候是冰凌灾害的重要基础信息，为冰凌洪水预测及河岸堤防修建提供重要依据［1］。河冰物候主要由四个关键物候节点来表征：1）初冰日FUS（freeze-up start）；2）完全冻结日FUE（freezeup end）；3）开始消融日BUS（break-up start）；4）完全消融日BUE（break-up end）［2-3］，这四个关键物候节点反映河流冰情的主要演变过程［4］。传统野外冰情观测需要较高的人力和物力成本，近些年遥感技术的多学科应用为冰凌物候研究及监测提供了全新手段［5-8］。当前，对于冰凌物候研究，多基于高时相分辨率的遥感数据，例如MODIS（Moderate-Resolution Imaging Spectroradiometer）及AVHRR（Advanced Very High Resolution Radiometer）等，且以海洋、湖泊等大尺度为主［2，9］。但这类数据较低空间分辨率限制了其在河道冰凌监测中的应用［10］。随着遥感影像时间和空间分辨率的不断提高，使得河冰表面物候提取成为可能。

目前，光学数据和雷达数据是两种主要用于监测大尺度冰凌物候信息的遥感数据源。对于光学数据，Nonaka等［9］基于MODIS数据对冰凌消融时间进行识别；Thuan等［3］基于MODIS产品利用反射率光谱曲线法识别了加拿大大奴湖物候信息，该方法能较为准确的反映河冰生消信息，但无法准确提取冰冻持续时间；Weber等［2］基于AVHRR数据开发了两步迭代的自动化算法，简化了湖冰物候信息提取过程；Li等［11］基于Landsat 7及MODIS等多源遥感数据，解译了阿蒙森海海冰过程及空间变异性。MODIS和AVHRR等大尺度光学数据在时间分辨率上有较大优势，保证了冰情监测的时效性，但其空间分辨率并不适用于黄河冰凌的监测。同时，Landsat及Sentinel系列等中高空间分辨率的光学数据在冬季易受云、雪覆盖影响，这使得单一传感器的光学数据的有效观测数量较少，往往难以满足物候监测的需求，而雷达数据不受雾、霾、烟尘等的阻挡，具有一定的穿透能力，提高了数据的可用性［12］。在雷达数据方面，Duguay等［13］基于RADARSAT-1监测湖冰生消相关过程，较为准确地识别了湖冰冻结与消融时间；Stephen等［14］基于QuickSCAT数据后向散射系数的时空演变，利用时序曲线法定义了湖冰的三个关键物候节点；Mermoz等［15］基于C波段RADARSAT-2数据开发了湖冰生消的三维模型，分析了湖冰生消的空间分布。这些研究多是以海洋、湖泊等大尺度为主的冰情信息监测，尚未应用于河冰物候的提取［16-17］。

当前，基于时间序列遥感数据提取冰凌物候信息的方法主要分为两大类，即阈值法和斜率法。阈值法大体分两种，第一种主要依据冰水间波谱特性提出的归一化差分积雪指数（NDSI，Normalized Difference Snow Index）提取河冰像元信息，根据冰期河冰覆盖率的变化人为划分阈值提取四个河冰物候节点［18-21］；这种方法在早期利用遥感技术监测冰情中较为常见，但要考虑河冰像元提取精度和阈值选取等问题。第二种主要根据凌汛期反射率光谱曲线或后向散射系数曲线划分阈值提取初冰日和完全消融日［22］；这种方法对于提取初冰日和完全消融日的提取精度较高，但仅适用于单一研究区或冰水后向散射系数相似的区域，对于范围跨度大、地形地貌差异性强的黄河存在难以拟定阈值的不足。类似的问题也出现在植被物候提取中，Jonsson等［23］2002年提出的动态阈值法对于上述问题的解决提供了参考。斜率法主要根据反射率时序曲线或后向散射系数曲线，将时序曲线的两个拐点分别定义为河冰的完全冻结日和开始消融日［24］。这种方法仅需要一组控制点的坐标信息就可以应用于多时相的遥感数据中，避免了人工提取中存在的偶发性误差和阈值法中因误提河冰像元而导致关键物候节点提取精度问题；考虑到时序曲线的完整性，该方法仅适用于提供逐日数据的遥感影像；为克服这个不足，Weber等［2］、Jeffries等［25］、Zhai等［26-27］基于曲线拟合的思想对凌汛期反射率光谱曲线或后向散射系数曲线进行重建，克服了数据缺失或不足的问题，在其研究区得到了较好的验证。综合来看，单一使用上述方法存在难以全面提取物候节点的问题，故两种方法结合使用有助于提高河冰物候提取精度。

凌汛灾害在我国北方以黄河（内蒙古段）最为典型，其孕灾环境复杂、致灾因子多样、链发性强，从而难于预测、防护以及抢险［28-29］。鉴于此，本文以黄河（内蒙古段）5个子段为研究对象，根据2015—2020年凌汛期时序Sentinel-1雷达数据的后向散射系数，利用斜率法及动态阈值法来提取5个子段河冰关键物候节点，以期为遥感技术监测河流冰情及研判冰凌洪水提供参考。

1 研究区及数据

1.1 研究区概况

研究区始于黄河内蒙古段海勃湾水库，终于万家寨水库库尾区，全长约为650 km，整体呈“几”字形分布［30］。河段冰期持续约100多天，封冻历时约20～30天，具有明显的物候差异，多年平均初冰日、完全冻结日及完全消融日分别约为11月22日、12月16日以及3月17日［31］。按照黄河防凌部署，开河期会存在人工破冰措施，为保证所研究的冰凌物候信息不受人工干扰，故选取的5个典型试验河段，分别为：海勃湾库尾区（S1）、三盛公水利枢纽闸下20 km内巴彦高勒河段（S2）、独贵塔拉镇三湖河口河段（S3）、什四份子弯道（S4）以及万家寨水库库尾河段（S5），5个试验段分别为库区型、过渡型、游荡型、弯曲型及峡谷型河段（图1），其中，S1段平均宽度为3 500 m，S2～S5段河宽介于200～600 m之间。库区型河段用于蓄水发电，深水区较多；过渡型河段浅滩多且河道摆动幅度大；游荡型河段淤积严重且经常发生河冰漫滩、偎堤现象；弯曲型河段易发生卡冰结坝现象；峡谷型则是黄河流出内蒙古段的特有河段［32］。

图1 研究区概况图Fig.1 Overview maps of the study area

1.2 数据来源与验证

1.2.1 Sentinel-1时间序列数据

Sentinel-1数据是欧洲航天局哥白尼计划（GMES）中的地球观测卫星，于2014年4月3日发射，载有C波段合成孔径雷达，可提供连续图像（白天、夜晚和各种天气），共具4种成像模式，分辨率最高5 m，幅宽达到400 km。本研究收集2015—2020年凌汛期所有覆盖试验区的Sentinel-1雷达影像，相对轨道号包含11、113、62、84、164以及157（图2），其中，S1～S5段2019—2020年54景、2018—2019年46景、2017—2018年47景、2016—2017年39景、2015—2016年21景，共计207景，源文件大小约5300GB。由于Sentinel-1数据2015年才逐步对外开放，故2015—2016年冰期数据量较少，2015—2016年除S1段外，所选取的数据已基本涵盖四个物候期信息。2019—2020年Sentinel-1数据来源于欧洲航天局（European Space Agency，https：//scihub.copernicus.eu/dhus/#home），由于历史数据在欧空局需逐个申请下载，故2015—2019年Sentinel-1数据从ASF DAAC数据库获得（https：//search.asf.alaska.edu/#/）。数据均为L1级IW模式下单视复数影像。由于近5年物候期和封冻历时等的不同，且一个凌汛期横跨两个年份，不利于物候期的定量描述，在此定义一种凌汛期内用于连续计算日期的方法：凌汛期累积日，即当年的11月10日至翌年3月31日（11月10日为1 d，以此类推）。

图2 Sentinel-1雷达数据Fig.2 radar data from Sentinel-1

1.2.2 验证数据

为保证物候信息提取的准确性及可靠性，本研究共使用两个数据源对物候信息进行验证：1）黄河网2015—2020年防凌专题（http：//www.yrcc.gov.cn/）凌汛期逐日观测数据；2）欧空局（https：//qc.sentinel1.eo.esa.int/aux_pp2/）2015—2020年所有覆盖研究区且云覆盖小于10%的Sentinel-2高分辨率光学遥感影像数据，共计264景，源文件大小约290 GB。其中，黄河网提供黄河（内蒙古段）各断面凌汛期流凌密度、封河位置及开河情况等书面信息（以S4段为例，表1），Sentinel-2高清影像作为视觉上物候信息的验证。Sentinel-2影像全部经过欧洲航天局提供的Sen2Cor（http：//step.esa.int/main/third-party-plugins-2/sen2cor/）插件处理成L2A级别大气底层反射率影像便于观测（以S1和S4段为例，图3）。验证数据以高分辨率光学影像为主，但当影像无法准确判断日期时，辅以黄河网凌汛期观测的书面数据。最终得到所有河段全部物候期作为验证数据。

表1 2015—2020年凌汛期S4段黄河网数据Table 1 Data of S4 section of the Yellow River Network during the ice flood season from 2015 to 2020

图3 2019—2020年S1和S5段河冰物候信息验证Fig.3 Verification of river ice phenological information in S1 and S5 sections from 2019 to 2020

2 研究方法

2.1 数据预处理及后向散射系数提取

Sentinel-1雷达数据按辐射定标、地理编码、多时相滤波、多时相配准以及裁剪进行预处理。由于数据量较大，单一时相滤波操作工作量过于繁冗，故选用多时相滤波同时操作多组数据，该操作在滤波的同时可将多时相的噪声进行统计，减少了后期统计分析的误差。辐射定标及地理编码在SNAP（Sentinel Application Platform）软 件（http：//step.esa.int/main/download）实现，该软件可实现两个步骤的归一化、批量化处理；导入雷达数据所需的精密轨道文件来自欧空局（https：//qc.sentinel1.eo.esa.int/aux_pp2/）；多时相配准所需的DEM数据为空间分辨率30 m的SRTM数据，来源于美国航空航天局（National Aeronautics and Space Administration简称NASA，http：//srtm.csi.cgiar.org/srtmdata/）；其余操作在SAR Scape工具箱中实现，最终目的是将SAR原始数据的DN（digital number）值转换为含有真实物理意义的后向散射系数值（σ0）用于物候分析。Sentinel-1雷达数据中共有两种极化方式：即VV和VH极化方式，Helean等［33］研究表明，VV极化下冰水阈值的差异性较VH极化下显著，故选用VV极化进行河冰物候研究。

后向散射系数表示微波能量返回功率与传输功率之比，一般将其取对数，用dB值表示，即：

且与地物的介电性能密切相关［34］。地物的介电性能描述为：

式中：ε'是介电常数；ε''为介电损耗，j=(-1)0.5。能量穿透深度（δp）可表示为：

式中：λ是波长；δp表示地物对σ0贡献的最大深度（如水体表现为强吸收，而冰雪则表现为强反射）。

S1～S5段多时相后向散射系数主要通过5组控制点进行提取。将预处理好的雷达影像导入SAR scape工具箱中，在各河段选取若干控制点。因各河段差异性较为明显，故控制点的数量并不做统一要求，但在选择控制点时应注意以下事项：1）控制点应均匀布设在各河段中，以便减少提取后向散射系数中的偶发性误差；2）由于清沟（在流速较大的急流、浅滩处或排水口处常出现狭长型不结冰的裸露水面）的存在，控制点应用各段稳封期影像选取，这样可减少后向散射系数提取的系统性误差。

2.2 时间序列数据拟合

河冰过程大体分为三个阶段，即生长阶段、稳定阶段以及消融阶段，σ0会伴随着这三个过程变化，从而使得物候曲线从生长过程到稳定过程呈“S”型，稳定过程到消融过程呈倒“S”型。常用分段式Logistic函数拟合、Savizky-Golay滤波器拟合以及Lowess平滑的方法［24］，由于S-G滤波要求时序数据连续且等间距，经过试验对比，分段Logistic方法拟合效果较另外两种拟合方法更能反映河冰物候规律，故选用分段Logistic拟合的方法来建立物候曲线。Logistic函数可以表示为：

式中：t为凌汛期累积日；f(t)为第t个凌汛期累积日的后向散射系数值；a和b为函数待求的拟合参数；d为河冰指数初始背景值，即：水体后向散射系数值；c+d为冰期冰凌指数最大值。

河冰生长阶段，σ0上升至某一点后在一定时期基本保持不变，以该范围的中点为分段依据进行分段。拟合的Logistic曲线共计48条，曲线经过若干次迭代后均收敛，误差趋近于最小，拟合优度均在0.9以上。

2.3 物候特征提取

河冰的σ0值在不同时期内发生的变化不同（变化速率和变化幅度等）［33］。生长阶段σ0值存在明显突变，变化速率快、变化幅度较大，不同河段变化幅度不同，无法统一阈值进行FUS和BUE的划分；稳定阶段σ0值增长速率放缓，FUE和BUS的阈值难以拟定［35］。动态阈值法可根据不同河段σ0变化幅度不同进行阈值划分，而斜率法则可在σ0增速放缓阶段找到其突变点。同时，考虑到传感器之间存在的差异（入射角和方位角等），因此本研究综合斜率法和动态阈值法提取黄河河冰物候节点。

2.3.1 FUE和BUS提取

河冰生长阶段向稳定阶段转换时，σ0增幅开始转变，Logistic函数曲线从一个近线性趋势转向另外一个近线性趋势，通过计算后向散射系数曲线的拐点，即极值可以确定河冰物候节点［2，24］，曲线曲率的公式可以根据公式（4）求导所得：

式中：z=exp(a+bt)；dα为曲线移动单位弧长时切线所转角度；ds为单位弧长。

结合河冰物候特性，Logistic函数曲线斜率最大值（最小值）为河冰生长（消融）速率变化最大时期，如图4（a）所示（以2019—2020年S2段为例）。本研究Logistic函数拟合的河冰生消曲线得到一个曲率最大值点［Logistic函数上升曲线，图4（a）］和一个曲率最小值点［Logistic函数下降曲线，图4（b）］，分别对应河冰物候期中FUE和BUS。

2.3.2 FUS和BUE提取

由于含沙量及光照等的差异性，黄河各段河冰和水体的后向散射系数值不尽相同。分别从各河段选取1 000个冰水像元进行分析（图5）。利用单一冰水间后向散射系数阈值差异性无法同时提取5段FUS和BUE。为了便于定性及定量化研究黄河河冰物候，故引入植被物候学中的动态阈值法（Sa）提取FUS和BUE［23］。动态阈值法指从基础水平开始计算，当拟合曲线的前半部分达到振幅指定的分数值，即定义为FUS；BUE由拟合曲线的后半部分达到指定的分数值来定义。数值选取不宜过大，否则会超过曲率最大值点与实际不符，同时也不宜过小，否则会使物候期提取误差过大。本文以0.1阈值为中心，均匀地设置了0.05、0.1、0.15及0.2四个阈值，简称为0.05 Sa、0.1 Sa、0.15 Sa及0.2 Sa（如，0.05 Sa表示从拟合曲线的下限开始，幅值达到曲线振幅的0.05记为0.05 Sa），如图4（c）所示。

图4 两种方法提取河冰关键物候节点Fig.4 Two methods to extract the key phenological nodes of river ice

图5 不同河段冰水像元后向散射系数值分布Fig.5 Numerical distribution of backscattering coefficients of ice water pixels in different river sections

2.4 精度评价

以2019—2020年凌汛期S1和S4段为例进行描述（图3）。根据黄河网逐日观测数据描述，S1段河冰四个关键物候节点分别为14 d、27 d、124 d及131 d，S4段分别为44 d、55 d、83 d以及106 d，通过光学遥感影像观测到的物候信息大体与该数据保持一致。同时，为便于描述河冰物候信息提取结果，采用偏差定量评价指标进行分析：

式中：RSi表示第i年河冰物候；Obseri表示验证数据在第i年的河冰物候。同时结合理论分析，定性描述试验段河冰物候特征及其变化规律。

3 结果与分析

3.1 河冰物候提取结果

表2总结了Sentinel-1后向散射系数2015—2020年在S1～S5段提取到的河冰物候变量。总体来看，基于动态阈值法提取到FUS和BUE的最大BIAS为+12 d、最小为0 d，而斜率法提取到BUE和FUS的BIAS最大为+5 d，最小为0 d，总体BIAS幅度为23 d。FUS从0.05 Sa～0.2 Sa的发展模式普遍为“延迟-提前”，而BUE的发展模式则为“提前-延迟”；FUS和BUE最优提取指标基本出现在“延迟-提前”之间的过渡阶段。同时发现，FUS和BUE对于阈值变化较为敏感，从0.05 Sa～0.1 Sa就能引起最小1 d、最大6 d的BIAS。FUE和BUS的提取结果规律性较弱，在48个提取结果中，BIAS为“提前”的为25个，占比52.1%，“延迟”占比47.9%。近5年FUE和BUS的提取精度逐渐增加，如S4段FUE的BIAS从+5 d逐渐提高至+1 d，而FUS的BIAS从-3 d逐渐提高至+2 d，可能与逐渐增多的数据有关。

表2 关键物候节点两种方法提取结果（BIAS）Table 2 Two methods to extract the results of key phenological nodes（BIAS）

3.1.1 FUE和BUS提取结果

由于各年份所获取的数据并非同一时相，且不均衡（图2），故需对不同数据量下提取到的河冰物候期进行分析。以S4段2015—2016年、2017—2018年及2019—2020年凌汛期拟合结果为例进行说明（图6）。基于曲线斜率法提取到不同数据量下的FUE和BUS，从拟合结果来看，不同数据量下拟合优度均在0.9以上，拟合效果较好。具体来看，FUS至FUE段（Logistic曲线前段），2015—2016年、2017—2018年、2019—2020年凌汛期拟合优度分别是0.96、0.92及0.94，2015年数据量最少，但其拟合优度最高，可能与样本的类型和特征有关［36］；BUS至BUE段（Logistic曲线后段）拟合优度分别为0.94、0.96及0.98，表明使用多点后向散射系数的均值进行曲线拟合比较符合河冰生长规律。从提取结果来看，2015—2016年FUE和BUS的BIAS分别为+5 d和-3 d，2017—2018年分别为-4 d和-2 d，而2019—2020年则是+1 d和+3 d，提取精度呈由低到高的趋势。同时，河冰的FUE和BUS并不是突发性事件，需要野外踏勘人工判别，因此，最大BIAS为+5 d精度良好［14］。

图6 S4段2015—2016年、2017—2018年及2019—2020年拟合结果Fig.6 Fitting results of S4 segment from 2015 to 2016，2017 to 2018 and 2019 to 2020

3.1.2 FUS和BUE提取结果

基于动态阈值法（Sa）提取到S1～S5段FUS和BUE，各河段的最优提取指标并不统一，但最优提取指标分别为0.1 Sa、0.2 Sa、0.1 Sa、0.05 Sa以及0.05 Sa，识 别 结 果 与Obser资 料 相 似，BIAS均在±3 d以内。同时，与2017—2018年和2019—2020年在其最优阈值提取到的河冰物候结果相比，2015—2016年所提取的精度较后两年数据差，其FUS提取精度较后两年数据充足情况下偏早2 d，BUE则偏晚2 d，可见物候节点提取精度与数据拟合结果和数据量有关。

3.2 河冰物候特征变化规律

3.2.1 河冰物候时间变化规律

黄河（内蒙古段）5个子段河冰物候节点均集中在10～128 d（11月中下旬至翌年2月中下旬，图7）S3～S4段FUS出 现 最 早，平 均 为14.7 d（11月24日），S1段出现最晚，平均约为38.5 d（12月18日）。2015—2020年各段FUS均有不同程度的偏晚现象，S1和S2段表现较为明显，分别偏晚约7 d和5 d；S1～S5段近5年FUS平均变化速率分别为1.4、1.0、0.8、0.2 d·a-1及0.4 d·a-1。S2段近5年FUS波动范围较大，波动幅度约15 d，其余各段波动幅度均在10 d左右。各段FUE主要集中在12月中旬至1月上旬，FUE最早和最晚分别是S3和S1段，近5年平均值分别约为24.6 d（11月24日）和48.2 d（12月28日）；除S2段，FUE波动范围均超过20 d，其余各段均在10 d左右。各段从FUS至FUE所用时间不尽相同，S1～S5段近5年平均约为9.7、14.0、8.2、12.8 d及10.8 d，其中S2段2016—2017年用时最长，为36 d。

图7 试验段关键物候节点时间分布Fig.7 Time distribution diagram of key phenological nodes in the test section

相比较而言，S1～S5段BUS较为分散，可从83 d（S1段）开始一直持续到126 d（S4段），大多子段的BUS相继出现在110～125 d，各段之间差异性较大，S1与S4段之间较为明显，两段BUS相差约22 d；除S4段，其余各段BUS均有不同偏早现象，变化速率分别为0.8、2.6、2.6 d·a-1及7.0 d·a-1。S1段近5年BUS波动范围较大，约27 d；S4段波动范围稳定在5 d以内，其余三段波动范围在15 d左右。各段BUE主要集中在三月，其中S1段最早（113 d），S4段最晚（129 d）。与BUS变化规律相似，除S1和S4段，其余三段BUE均有不同程度偏早现象，变化速率分别为1.6、1.0 d·a-1及7.2 d·a-1。在研究时段内，5段由BUE至FUS的平均封冻期为77 d，其中S1段封冻期最短（47 d），S4段封冻期最长（96 d）。

进一步分析S1～S5段物候信息发现，S1和S4段从BUS至BUE的 横向冻 结速 率 正以0.2 d·a-1和1.4 d·a-1的速度逐年增加，而S3和S5段的横向冻结速率呈1.4 d·a-1和1.0 d·a-1的速度逐年减少，S2段近5年基本保持不变，由于每年FUS均集中发生于S3～S5段，表明黄河内蒙古段近5年首封位置逐步向下游偏移。S1～S3段的横向消融速率呈增加趋势，分别为0.4、1.0、1.8 d·a-1，S4段基本保持不变，而S5段以0.6 d·a-1横向消融速度呈放缓趋势。

3.2.2 河冰物候冻结（消融）空间模式

图8～9给出了2015—2020年试验段物候期FUS及BUE的冻结（消融）空间模式。总体来讲，研究时间内5段间FUS和BUE存在较大差异；对于FUS，一般S1和S2段最先流凌，然后逐步向其他各段延伸，总体呈“中间早，两边晚”的发展模式；FUE的发展模式与FUS相反，大体呈“两边早，中间晚”的格局，表明黄河（内蒙古段）封开河逆序且各段逐步发展的空间模式。以S1～S5段进行分析，S1段河冰FUS发展模式为“两岸浅水区—库区中心—库尾深水区”，而BUE则呈“库区中心—两岸浅水区及深水区”的模式发展；S2～S5段相似，发展模式为“两岸或河滩—主流区”；需要说明的是，弯道处FUS的发展模式为“凹岸—凸岸—主流区”，而BUE则按照“主流区—凸岸—凹岸”的顺序逐步发展。

图8 试验段FUS时空分布Fig.8 Temporal and spatial distribution of FUS in test section

对S1～S5段进一步分析发现，FUS和BUE的空间分布可以反映出水深及河底地形地貌的差异性。以S1段为例，张磊等［37］的研究表明，该河段的河底地势由东南向西北倾斜，深水区位于距坝体约860 m处的库尾段，浅水区位于两岸及西南部，结合S1段FUS可知，河冰的FUS空间模式与河底地形地貌分布具有一致性，而BUE的空间分布与FUS恰好相反。

4 讨论

冬季，由于河道上空雾、云及霾的存在，光学数据在识别河冰物候信息中受到限制，合成孔径雷达数据在物候识别方面效果较好［33］。同时，黄河河冰物候研究受限于数据的时相分辨率，因此，常在研究中出现某一物候期缺失数据的情况［18］。与大尺度海冰（湖冰）物候监测和植被物候监测不同，本文基于Sentinel-1时间序列数据，综合曲线斜率法和动态阈值法构建了黄河（内蒙古段）5个子段的河冰物候曲线，提出了S1～S5段河冰关键物候节点的最优提取指标，解决了河冰物候信息难以提取的问题［38］。本文发现，时序数据提取、拟合以及提取方法是影响提取精度的重要因素。在时间序列后向散射系数提取中，由于河道清沟和湿雪的存在，易将河冰误提取成水体，可能导致物候提取的BIAS较大，故应选用稳封期影像进行控制点的选取［25，34］。在数据拟合时，与S-G滤波、Lowess平滑等重构河冰物候曲线的方法相比，Logistic拟合更能反映河冰生消规律，同时不受时序数据等间距的限制，适用于黄河冰情解译［14，24］。各河段数据量不同，但拟合效果良好，主要与样本类型及特征有关［36，39］。在提取方法方面，由于各段后向散射系数的差异，单一阈值易造成FUS和BUE提取的失真，动态阈值根据各段物候曲线特性提取物候期，但无法量化导致各段之间差异性的不确定性因素（如含沙量、反射率、太阳高度角、河底地质构造等）。

黄河（内蒙古段）各子河段物候特征与前人湖冰物候特征的研究相比，既具统一性，也具差异性。S1段是典型库区型河段，水流流速相对较小，FUS和BUE的发展模式与湖冰发展模式基本相同［20］。S2～S5段顺直部分的FUS和BUE发展模式与S1段相似，区别在于，弯道处由于水流离心力的存在导致河冰生长多集中于凹岸侧，其FUS发展模式为“凹岸—凸岸—主流区”。同时，弯道凹岸处由于存在大量堆冰，使得河岸滩在冻融侵蚀和淘刷作用下易崩塌，故弯道凹岸侧应加强防护［40-42］。

黄河河冰物候特征不仅受气温、风速、太阳辐射以及积雪覆盖等的影响，并且也受水位、流量，河道边界条件以及上下游水库防凌调度的影响［43-45］，其中气温和水力条件被认为是影响黄河河冰物候信息的主要因子［46］。例如：气温越高初冰日日期越晚；水位越高，河道与大气表面接触面积就越大，河冰冻结所需要的热量条件也随之增大，从而导致河流消融时间变得缓慢等［47］。因此，在未来研究中，如何在模型中显示结合气象和水文数据是一个值得考虑的方向，这可能会进一步提高河冰表面物候的提取精度。

图9 试验段BUE时空分布Fig.9 Temporal and spatial distribution of BUE in test section

5 结论

本文基于2015—2020年间Sentinel-1/2雷达/光学影像数据，利用Logistic函数拟合，综合曲线斜率法及动态阈值法确定了5个试验段河冰关键物候节点的阈值，分析了河冰物候特征变化规律及其主要驱动因素，得到如下结论：

（1）基于Sentinel-1时间序列数据，综合斜率法和动态阈值法可用于黄河（内蒙古段）河冰物候提取，整体提取精度较高。

（2）S1～S5段FUS和BUE的河冰物候最优阈值提取指标分别为：0.1 Sa、0.2 Sa、0.1 Sa、0.05 Sa以及0.05 Sa，BIAS在±3 d以内；斜率法提取FUE和BUS的BIAS在±5 d以内。

（3）S1～S5段近5年FUS平均变化速率分别为1.4、1.0、0.8、0.2 d·a-1及0.4 d·a-1；S2段 近5年FUS波动范围较大，波动幅度为15 d左右，其余各段波动幅度均在10 d左右；除S4段，其余各段BUS均有不同程度偏早现象，变化速率分别为0.8、2.6、2.6 d·a-1及7 d·a-1；S1段近5年BUS波动范围较大，约为27 d；S4段波动范围稳定在5 d以内，其余三段波动范围在15 d左右。

（4）S1段FUS发展模式为“两岸浅水区—库区中心—库尾深水区”；S2～S5段FUS发展模式为“两岸或河滩—主流区”；弯道处FUS发展模式为“凹岸—凸岸—主流区”；S1段BUE发展模式为“库区中心—两岸浅水区及深水区”；S2～S5段顺直河段BUE发展模式为“主流区—两岸或河滩”；弯道处BUE发展模式为“主流区—凸岸—凹岸”。