基于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數據開發了湖冰生消的三維模型，分析了湖冰生消的空間分布。這些研究多是以海洋、湖泊等大尺度為主的冰情信息監測，尚未應用于河冰物候的提?。?6-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段提取到的河冰物候變量?？傮w來看，基于動態閾值法提取到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較大，故應選用穩封期影像進行控制點的選?。?5，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發展模式為“主流區—凸岸—凹岸”。