張 鑫，牟獻(xiàn)友，翟涌光，張寶森，冀鴻蘭

（1.內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院，內(nèi)蒙古 呼和浩特 010018；2.黃河水利委員會(huì)黃河水利科學(xué)研究院，河南 鄭州 450003）

0 引言

冰凌洪水災(zāi)害是高緯度河流特有的一種易發(fā)自然災(zāi)害，河冰物候是冰凌災(zāi)害的重要基礎(chǔ)信息，為冰凌洪水預(yù)測(cè)及河岸堤防修建提供重要依據(jù)［1］。河冰物候主要由四個(gè)關(guān)鍵物候節(jié)點(diǎn)來表征：1）初冰日FUS（freeze-up start）；2）完全凍結(jié)日FUE（freezeup end）；3）開始消融日BUS（break-up start）；4）完全消融日BUE（break-up end）［2-3］，這四個(gè)關(guān)鍵物候節(jié)點(diǎn)反映河流冰情的主要演變過程［4］。傳統(tǒng)野外冰情觀測(cè)需要較高的人力和物力成本，近些年遙感技術(shù)的多學(xué)科應(yīng)用為冰凌物候研究及監(jiān)測(cè)提供了全新手段［5-8］。當(dāng)前，對(duì)于冰凌物候研究，多基于高時(shí)相分辨率的遙感數(shù)據(jù)，例如MODIS（Moderate-Resolution Imaging Spectroradiometer）及AVHRR（Advanced Very High Resolution Radiometer）等，且以海洋、湖泊等大尺度為主［2，9］。但這類數(shù)據(jù)較低空間分辨率限制了其在河道冰凌監(jiān)測(cè)中的應(yīng)用［10］。隨著遙感影像時(shí)間和空間分辨率的不斷提高，使得河冰表面物候提取成為可能。

目前，光學(xué)數(shù)據(jù)和雷達(dá)數(shù)據(jù)是兩種主要用于監(jiān)測(cè)大尺度冰凌物候信息的遙感數(shù)據(jù)源。對(duì)于光學(xué)數(shù)據(jù)，Nonaka等［9］基于MODIS數(shù)據(jù)對(duì)冰凌消融時(shí)間進(jìn)行識(shí)別；Thuan等［3］基于MODIS產(chǎn)品利用反射率光譜曲線法識(shí)別了加拿大大奴湖物候信息，該方法能較為準(zhǔn)確的反映河冰生消信息，但無法準(zhǔn)確提取冰凍持續(xù)時(shí)間；Weber等［2］基于AVHRR數(shù)據(jù)開發(fā)了兩步迭代的自動(dòng)化算法，簡(jiǎn)化了湖冰物候信息提取過程；Li等［11］基于Landsat 7及MODIS等多源遙感數(shù)據(jù)，解譯了阿蒙森海海冰過程及空間變異性。MODIS和AVHRR等大尺度光學(xué)數(shù)據(jù)在時(shí)間分辨率上有較大優(yōu)勢(shì)，保證了冰情監(jiān)測(cè)的時(shí)效性，但其空間分辨率并不適用于黃河冰凌的監(jiān)測(cè)。同時(shí)，Landsat及Sentinel系列等中高空間分辨率的光學(xué)數(shù)據(jù)在冬季易受云、雪覆蓋影響，這使得單一傳感器的光學(xué)數(shù)據(jù)的有效觀測(cè)數(shù)量較少，往往難以滿足物候監(jiān)測(cè)的需求，而雷達(dá)數(shù)據(jù)不受霧、霾、煙塵等的阻擋，具有一定的穿透能力，提高了數(shù)據(jù)的可用性［12］。在雷達(dá)數(shù)據(jù)方面，Duguay等［13］基于RADARSAT-1監(jiān)測(cè)湖冰生消相關(guān)過程，較為準(zhǔn)確地識(shí)別了湖冰凍結(jié)與消融時(shí)間；Stephen等［14］基于QuickSCAT數(shù)據(jù)后向散射系數(shù)的時(shí)空演變，利用時(shí)序曲線法定義了湖冰的三個(gè)關(guān)鍵物候節(jié)點(diǎn)；Mermoz等［15］基于C波段RADARSAT-2數(shù)據(jù)開發(fā)了湖冰生消的三維模型，分析了湖冰生消的空間分布。這些研究多是以海洋、湖泊等大尺度為主的冰情信息監(jiān)測(cè)，尚未應(yīng)用于河冰物候的提取［16-17］。

當(dāng)前，基于時(shí)間序列遙感數(shù)據(jù)提取冰凌物候信息的方法主要分為兩大類，即閾值法和斜率法。閾值法大體分兩種，第一種主要依據(jù)冰水間波譜特性提出的歸一化差分積雪指數(shù)（NDSI，Normalized Difference Snow Index）提取河冰像元信息，根據(jù)冰期河冰覆蓋率的變化人為劃分閾值提取四個(gè)河冰物候節(jié)點(diǎn)［18-21］；這種方法在早期利用遙感技術(shù)監(jiān)測(cè)冰情中較為常見，但要考慮河冰像元提取精度和閾值選取等問題。第二種主要根據(jù)凌汛期反射率光譜曲線或后向散射系數(shù)曲線劃分閾值提取初冰日和完全消融日［22］；這種方法對(duì)于提取初冰日和完全消融日的提取精度較高，但僅適用于單一研究區(qū)或冰水后向散射系數(shù)相似的區(qū)域，對(duì)于范圍跨度大、地形地貌差異性強(qiáng)的黃河存在難以擬定閾值的不足。類似的問題也出現(xiàn)在植被物候提取中，Jonsson等［23］2002年提出的動(dòng)態(tài)閾值法對(duì)于上述問題的解決提供了參考。斜率法主要根據(jù)反射率時(shí)序曲線或后向散射系數(shù)曲線，將時(shí)序曲線的兩個(gè)拐點(diǎn)分別定義為河冰的完全凍結(jié)日和開始消融日［24］。這種方法僅需要一組控制點(diǎn)的坐標(biāo)信息就可以應(yīng)用于多時(shí)相的遙感數(shù)據(jù)中，避免了人工提取中存在的偶發(fā)性誤差和閾值法中因誤提河冰像元而導(dǎo)致關(guān)鍵物候節(jié)點(diǎn)提取精度問題；考慮到時(shí)序曲線的完整性，該方法僅適用于提供逐日數(shù)據(jù)的遙感影像；為克服這個(gè)不足，Weber等［2］、Jeffries等［25］、Zhai等［26-27］基于曲線擬合的思想對(duì)凌汛期反射率光譜曲線或后向散射系數(shù)曲線進(jìn)行重建，克服了數(shù)據(jù)缺失或不足的問題，在其研究區(qū)得到了較好的驗(yàn)證。綜合來看，單一使用上述方法存在難以全面提取物候節(jié)點(diǎn)的問題，故兩種方法結(jié)合使用有助于提高河冰物候提取精度。

凌汛災(zāi)害在我國(guó)北方以黃河（內(nèi)蒙古段）最為典型，其孕災(zāi)環(huán)境復(fù)雜、致災(zāi)因子多樣、鏈發(fā)性強(qiáng)，從而難于預(yù)測(cè)、防護(hù)以及搶險(xiǎn)［28-29］。鑒于此，本文以黃河（內(nèi)蒙古段）5個(gè)子段為研究對(duì)象，根據(jù)2015—2020年凌汛期時(shí)序Sentinel-1雷達(dá)數(shù)據(jù)的后向散射系數(shù)，利用斜率法及動(dòng)態(tài)閾值法來提取5個(gè)子段河冰關(guān)鍵物候節(jié)點(diǎn)，以期為遙感技術(shù)監(jiān)測(cè)河流冰情及研判冰凌洪水提供參考。

1 研究區(qū)及數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)始于黃河內(nèi)蒙古段海勃灣水庫，終于萬家寨水庫庫尾區(qū)，全長(zhǎng)約為650 km，整體呈“幾”字形分布［30］。河段冰期持續(xù)約100多天，封凍歷時(shí)約20～30天，具有明顯的物候差異，多年平均初冰日、完全凍結(jié)日及完全消融日分別約為11月22日、12月16日以及3月17日［31］。按照黃河防凌部署，開河期會(huì)存在人工破冰措施，為保證所研究的冰凌物候信息不受人工干擾，故選取的5個(gè)典型試驗(yàn)河段，分別為：海勃灣庫尾區(qū)（S1）、三盛公水利樞紐閘下20 km內(nèi)巴彥高勒河段（S2）、獨(dú)貴塔拉鎮(zhèn)三湖河口河段（S3）、什四份子彎道（S4）以及萬家寨水庫庫尾河段（S5），5個(gè)試驗(yàn)段分別為庫區(qū)型、過渡型、游蕩型、彎曲型及峽谷型河段（圖1），其中，S1段平均寬度為3 500 m，S2～S5段河寬介于200～600 m之間。庫區(qū)型河段用于蓄水發(fā)電，深水區(qū)較多；過渡型河段淺灘多且河道擺動(dòng)幅度大；游蕩型河段淤積嚴(yán)重且經(jīng)常發(fā)生河冰漫灘、偎堤現(xiàn)象；彎曲型河段易發(fā)生卡冰結(jié)壩現(xiàn)象；峽谷型則是黃河流出內(nèi)蒙古段的特有河段［32］。

圖1 研究區(qū)概況圖Fig.1 Overview maps of the study area

1.2 數(shù)據(jù)來源與驗(yàn)證

1.2.1 Sentinel-1時(shí)間序列數(shù)據(jù)

Sentinel-1數(shù)據(jù)是歐洲航天局哥白尼計(jì)劃（GMES）中的地球觀測(cè)衛(wèi)星，于2014年4月3日發(fā)射，載有C波段合成孔徑雷達(dá)，可提供連續(xù)圖像（白天、夜晚和各種天氣），共具4種成像模式，分辨率最高5 m，幅寬達(dá)到400 km。本研究收集2015—2020年凌汛期所有覆蓋試驗(yàn)區(qū)的Sentinel-1雷達(dá)影像，相對(duì)軌道號(hào)包含11、113、62、84、164以及157（圖2），其中，S1～S5段2019—2020年54景、2018—2019年46景、2017—2018年47景、2016—2017年39景、2015—2016年21景，共計(jì)207景，源文件大小約5300GB。由于Sentinel-1數(shù)據(jù)2015年才逐步對(duì)外開放，故2015—2016年冰期數(shù)據(jù)量較少，2015—2016年除S1段外，所選取的數(shù)據(jù)已基本涵蓋四個(gè)物候期信息。2019—2020年Sentinel-1數(shù)據(jù)來源于歐洲航天局（European Space Agency，https：//scihub.copernicus.eu/dhus/#home），由于歷史數(shù)據(jù)在歐空局需逐個(gè)申請(qǐng)下載，故2015—2019年Sentinel-1數(shù)據(jù)從ASF DAAC數(shù)據(jù)庫獲得（https：//search.asf.alaska.edu/#/）。數(shù)據(jù)均為L(zhǎng)1級(jí)IW模式下單視復(fù)數(shù)影像。由于近5年物候期和封凍歷時(shí)等的不同，且一個(gè)凌汛期橫跨兩個(gè)年份，不利于物候期的定量描述，在此定義一種凌汛期內(nèi)用于連續(xù)計(jì)算日期的方法：凌汛期累積日，即當(dāng)年的11月10日至翌年3月31日（11月10日為1 d，以此類推）。

圖2 Sentinel-1雷達(dá)數(shù)據(jù)Fig.2 radar data from Sentinel-1

1.2.2 驗(yàn)證數(shù)據(jù)

為保證物候信息提取的準(zhǔn)確性及可靠性，本研究共使用兩個(gè)數(shù)據(jù)源對(duì)物候信息進(jìn)行驗(yàn)證：1）黃河網(wǎng)2015—2020年防凌專題（http：//www.yrcc.gov.cn/）凌汛期逐日觀測(cè)數(shù)據(jù)；2）歐空局（https：//qc.sentinel1.eo.esa.int/aux_pp2/）2015—2020年所有覆蓋研究區(qū)且云覆蓋小于10%的Sentinel-2高分辨率光學(xué)遙感影像數(shù)據(jù)，共計(jì)264景，源文件大小約290 GB。其中，黃河網(wǎng)提供黃河（內(nèi)蒙古段）各斷面凌汛期流凌密度、封河位置及開河情況等書面信息（以S4段為例，表1），Sentinel-2高清影像作為視覺上物候信息的驗(yàn)證。Sentinel-2影像全部經(jīng)過歐洲航天局提供的Sen2Cor（http：//step.esa.int/main/third-party-plugins-2/sen2cor/）插件處理成L2A級(jí)別大氣底層反射率影像便于觀測(cè)（以S1和S4段為例，圖3）。驗(yàn)證數(shù)據(jù)以高分辨率光學(xué)影像為主，但當(dāng)影像無法準(zhǔn)確判斷日期時(shí)，輔以黃河網(wǎng)凌汛期觀測(cè)的書面數(shù)據(jù)。最終得到所有河段全部物候期作為驗(yàn)證數(shù)據(jù)。

表1 2015—2020年凌汛期S4段黃河網(wǎng)數(shù)據(jù)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段河冰物候信息驗(yàn)證Fig.3 Verification of river ice phenological information in S1 and S5 sections from 2019 to 2020

2 研究方法

2.1 數(shù)據(jù)預(yù)處理及后向散射系數(shù)提取

Sentinel-1雷達(dá)數(shù)據(jù)按輻射定標(biāo)、地理編碼、多時(shí)相濾波、多時(shí)相配準(zhǔn)以及裁剪進(jìn)行預(yù)處理。由于數(shù)據(jù)量較大，單一時(shí)相濾波操作工作量過于繁冗，故選用多時(shí)相濾波同時(shí)操作多組數(shù)據(jù)，該操作在濾波的同時(shí)可將多時(shí)相的噪聲進(jìn)行統(tǒng)計(jì)，減少了后期統(tǒng)計(jì)分析的誤差。輻射定標(biāo)及地理編碼在SNAP（Sentinel Application Platform）軟 件（http：//step.esa.int/main/download）實(shí)現(xiàn)，該軟件可實(shí)現(xiàn)兩個(gè)步驟的歸一化、批量化處理；導(dǎo)入雷達(dá)數(shù)據(jù)所需的精密軌道文件來自歐空局（https：//qc.sentinel1.eo.esa.int/aux_pp2/）；多時(shí)相配準(zhǔn)所需的DEM數(shù)據(jù)為空間分辨率30 m的SRTM數(shù)據(jù)，來源于美國(guó)航空航天局（National Aeronautics and Space Administration簡(jiǎn)稱NASA，http：//srtm.csi.cgiar.org/srtmdata/）；其余操作在SAR Scape工具箱中實(shí)現(xiàn)，最終目的是將SAR原始數(shù)據(jù)的DN（digital number）值轉(zhuǎn)換為含有真實(shí)物理意義的后向散射系數(shù)值（σ0）用于物候分析。Sentinel-1雷達(dá)數(shù)據(jù)中共有兩種極化方式：即VV和VH極化方式，Helean等［33］研究表明，VV極化下冰水閾值的差異性較VH極化下顯著，故選用VV極化進(jìn)行河冰物候研究。

后向散射系數(shù)表示微波能量返回功率與傳輸功率之比，一般將其取對(duì)數(shù)，用dB值表示，即：

且與地物的介電性能密切相關(guān)［34］。地物的介電性能描述為：

式中：ε'是介電常數(shù)；ε''為介電損耗，j=(-1)0.5。能量穿透深度（δp）可表示為：

式中：λ是波長(zhǎng)；δp表示地物對(duì)σ0貢獻(xiàn)的最大深度（如水體表現(xiàn)為強(qiáng)吸收，而冰雪則表現(xiàn)為強(qiáng)反射）。

S1～S5段多時(shí)相后向散射系數(shù)主要通過5組控制點(diǎn)進(jìn)行提取。將預(yù)處理好的雷達(dá)影像導(dǎo)入SAR scape工具箱中，在各河段選取若干控制點(diǎn)。因各河段差異性較為明顯，故控制點(diǎn)的數(shù)量并不做統(tǒng)一要求，但在選擇控制點(diǎn)時(shí)應(yīng)注意以下事項(xiàng)：1）控制點(diǎn)應(yīng)均勻布設(shè)在各河段中，以便減少提取后向散射系數(shù)中的偶發(fā)性誤差；2）由于清溝（在流速較大的急流、淺灘處或排水口處常出現(xiàn)狹長(zhǎng)型不結(jié)冰的裸露水面）的存在，控制點(diǎn)應(yīng)用各段穩(wěn)封期影像選取，這樣可減少后向散射系數(shù)提取的系統(tǒng)性誤差。

2.2 時(shí)間序列數(shù)據(jù)擬合

河冰過程大體分為三個(gè)階段，即生長(zhǎng)階段、穩(wěn)定階段以及消融階段，σ0會(huì)伴隨著這三個(gè)過程變化，從而使得物候曲線從生長(zhǎng)過程到穩(wěn)定過程呈“S”型，穩(wěn)定過程到消融過程呈倒“S”型。常用分段式Logistic函數(shù)擬合、Savizky-Golay濾波器擬合以及Lowess平滑的方法［24］，由于S-G濾波要求時(shí)序數(shù)據(jù)連續(xù)且等間距，經(jīng)過試驗(yàn)對(duì)比，分段Logistic方法擬合效果較另外兩種擬合方法更能反映河冰物候規(guī)律，故選用分段Logistic擬合的方法來建立物候曲線。Logistic函數(shù)可以表示為：

式中：t為凌汛期累積日；f(t)為第t個(gè)凌汛期累積日的后向散射系數(shù)值；a和b為函數(shù)待求的擬合參數(shù)；d為河冰指數(shù)初始背景值，即：水體后向散射系數(shù)值；c+d為冰期冰凌指數(shù)最大值。

河冰生長(zhǎng)階段，σ0上升至某一點(diǎn)后在一定時(shí)期基本保持不變，以該范圍的中點(diǎn)為分段依據(jù)進(jìn)行分段。擬合的Logistic曲線共計(jì)48條，曲線經(jīng)過若干次迭代后均收斂，誤差趨近于最小，擬合優(yōu)度均在0.9以上。

2.3 物候特征提取

河冰的σ0值在不同時(shí)期內(nèi)發(fā)生的變化不同（變化速率和變化幅度等）［33］。生長(zhǎng)階段σ0值存在明顯突變，變化速率快、變化幅度較大，不同河段變化幅度不同，無法統(tǒng)一閾值進(jìn)行FUS和BUE的劃分；穩(wěn)定階段σ0值增長(zhǎng)速率放緩，F(xiàn)UE和BUS的閾值難以擬定［35］。動(dòng)態(tài)閾值法可根據(jù)不同河段σ0變化幅度不同進(jìn)行閾值劃分，而斜率法則可在σ0增速放緩階段找到其突變點(diǎn)。同時(shí)，考慮到傳感器之間存在的差異（入射角和方位角等），因此本研究綜合斜率法和動(dòng)態(tài)閾值法提取黃河河冰物候節(jié)點(diǎn)。

2.3.1 FUE和BUS提取

河冰生長(zhǎng)階段向穩(wěn)定階段轉(zhuǎn)換時(shí)，σ0增幅開始轉(zhuǎn)變，Logistic函數(shù)曲線從一個(gè)近線性趨勢(shì)轉(zhuǎn)向另外一個(gè)近線性趨勢(shì)，通過計(jì)算后向散射系數(shù)曲線的拐點(diǎn)，即極值可以確定河冰物候節(jié)點(diǎn)［2，24］，曲線曲率的公式可以根據(jù)公式（4）求導(dǎo)所得：

式中：z=exp(a+bt)；dα為曲線移動(dòng)單位弧長(zhǎng)時(shí)切線所轉(zhuǎn)角度；ds為單位弧長(zhǎng)。

結(jié)合河冰物候特性，Logistic函數(shù)曲線斜率最大值（最小值）為河冰生長(zhǎng)（消融）速率變化最大時(shí)期，如圖4（a）所示（以2019—2020年S2段為例）。本研究Logistic函數(shù)擬合的河冰生消曲線得到一個(gè)曲率最大值點(diǎn)［Logistic函數(shù)上升曲線，圖4（a）］和一個(gè)曲率最小值點(diǎn)［Logistic函數(shù)下降曲線，圖4（b）］，分別對(duì)應(yīng)河冰物候期中FUE和BUS。

2.3.2 FUS和BUE提取

由于含沙量及光照等的差異性，黃河各段河冰和水體的后向散射系數(shù)值不盡相同。分別從各河段選取1 000個(gè)冰水像元進(jìn)行分析（圖5）。利用單一冰水間后向散射系數(shù)閾值差異性無法同時(shí)提取5段FUS和BUE。為了便于定性及定量化研究黃河河冰物候，故引入植被物候?qū)W中的動(dòng)態(tài)閾值法（Sa）提取FUS和BUE［23］。動(dòng)態(tài)閾值法指從基礎(chǔ)水平開始計(jì)算，當(dāng)擬合曲線的前半部分達(dá)到振幅指定的分?jǐn)?shù)值，即定義為FUS；BUE由擬合曲線的后半部分達(dá)到指定的分?jǐn)?shù)值來定義。數(shù)值選取不宜過大，否則會(huì)超過曲率最大值點(diǎn)與實(shí)際不符，同時(shí)也不宜過小，否則會(huì)使物候期提取誤差過大。本文以0.1閾值為中心，均勻地設(shè)置了0.05、0.1、0.15及0.2四個(gè)閾值，簡(jiǎn)稱為0.05 Sa、0.1 Sa、0.15 Sa及0.2 Sa（如，0.05 Sa表示從擬合曲線的下限開始，幅值達(dá)到曲線振幅的0.05記為0.05 Sa），如圖4（c）所示。

圖4 兩種方法提取河冰關(guān)鍵物候節(jié)點(diǎn)Fig.4 Two methods to extract the key phenological nodes of river ice

圖5 不同河段冰水像元后向散射系數(shù)值分布Fig.5 Numerical distribution of backscattering coefficients of ice water pixels in different river sections

2.4 精度評(píng)價(jià)

以2019—2020年凌汛期S1和S4段為例進(jìn)行描述（圖3）。根據(jù)黃河網(wǎng)逐日觀測(cè)數(shù)據(jù)描述，S1段河冰四個(gè)關(guān)鍵物候節(jié)點(diǎn)分別為14 d、27 d、124 d及131 d，S4段分別為44 d、55 d、83 d以及106 d，通過光學(xué)遙感影像觀測(cè)到的物候信息大體與該數(shù)據(jù)保持一致。同時(shí)，為便于描述河冰物候信息提取結(jié)果，采用偏差定量評(píng)價(jià)指標(biāo)進(jìn)行分析：

式中：RSi表示第i年河冰物候；Obseri表示驗(yàn)證數(shù)據(jù)在第i年的河冰物候。同時(shí)結(jié)合理論分析，定性描述試驗(yàn)段河冰物候特征及其變化規(guī)律。

3 結(jié)果與分析

3.1 河冰物候提取結(jié)果

表2總結(jié)了Sentinel-1后向散射系數(shù)2015—2020年在S1～S5段提取到的河冰物候變量。總體來看，基于動(dòng)態(tài)閾值法提取到FUS和BUE的最大BIAS為+12 d、最小為0 d，而斜率法提取到BUE和FUS的BIAS最大為+5 d，最小為0 d，總體BIAS幅度為23 d。FUS從0.05 Sa～0.2 Sa的發(fā)展模式普遍為“延遲-提前”，而BUE的發(fā)展模式則為“提前-延遲”；FUS和BUE最優(yōu)提取指標(biāo)基本出現(xiàn)在“延遲-提前”之間的過渡階段。同時(shí)發(fā)現(xiàn)，F(xiàn)US和BUE對(duì)于閾值變化較為敏感，從0.05 Sa～0.1 Sa就能引起最小1 d、最大6 d的BIAS。FUE和BUS的提取結(jié)果規(guī)律性較弱，在48個(gè)提取結(jié)果中，BIAS為“提前”的為25個(gè)，占比52.1%，“延遲”占比47.9%。近5年FUE和BUS的提取精度逐漸增加，如S4段FUE的BIAS從+5 d逐漸提高至+1 d，而FUS的BIAS從-3 d逐漸提高至+2 d，可能與逐漸增多的數(shù)據(jù)有關(guān)。

表2 關(guān)鍵物候節(jié)點(diǎn)兩種方法提取結(jié)果（BIAS）Table 2 Two methods to extract the results of key phenological nodes（BIAS）

3.1.1 FUE和BUS提取結(jié)果

由于各年份所獲取的數(shù)據(jù)并非同一時(shí)相，且不均衡（圖2），故需對(duì)不同數(shù)據(jù)量下提取到的河冰物候期進(jìn)行分析。以S4段2015—2016年、2017—2018年及2019—2020年凌汛期擬合結(jié)果為例進(jìn)行說明（圖6）。基于曲線斜率法提取到不同數(shù)據(jù)量下的FUE和BUS，從擬合結(jié)果來看，不同數(shù)據(jù)量下擬合優(yōu)度均在0.9以上，擬合效果較好。具體來看，F(xiàn)US至FUE段（Logistic曲線前段），2015—2016年、2017—2018年、2019—2020年凌汛期擬合優(yōu)度分別是0.96、0.92及0.94，2015年數(shù)據(jù)量最少，但其擬合優(yōu)度最高，可能與樣本的類型和特征有關(guān)［36］；BUS至BUE段（Logistic曲線后段）擬合優(yōu)度分別為0.94、0.96及0.98，表明使用多點(diǎn)后向散射系數(shù)的均值進(jìn)行曲線擬合比較符合河冰生長(zhǎng)規(guī)律。從提取結(jié)果來看，2015—2016年FUE和BUS的BIAS分別為+5 d和-3 d，2017—2018年分別為-4 d和-2 d，而2019—2020年則是+1 d和+3 d，提取精度呈由低到高的趨勢(shì)。同時(shí)，河冰的FUE和BUS并不是突發(fā)性事件，需要野外踏勘人工判別，因此，最大BIAS為+5 d精度良好［14］。

圖6 S4段2015—2016年、2017—2018年及2019—2020年擬合結(jié)果Fig.6 Fitting results of S4 segment from 2015 to 2016，2017 to 2018 and 2019 to 2020

3.1.2 FUS和BUE提取結(jié)果

基于動(dòng)態(tài)閾值法（Sa）提取到S1～S5段FUS和BUE，各河段的最優(yōu)提取指標(biāo)并不統(tǒng)一，但最優(yōu)提取指標(biāo)分別為0.1 Sa、0.2 Sa、0.1 Sa、0.05 Sa以及0.05 Sa，識(shí) 別 結(jié) 果 與Obser資 料 相 似，BIAS均在±3 d以內(nèi)。同時(shí)，與2017—2018年和2019—2020年在其最優(yōu)閾值提取到的河冰物候結(jié)果相比，2015—2016年所提取的精度較后兩年數(shù)據(jù)差，其FUS提取精度較后兩年數(shù)據(jù)充足情況下偏早2 d，BUE則偏晚2 d，可見物候節(jié)點(diǎn)提取精度與數(shù)據(jù)擬合結(jié)果和數(shù)據(jù)量有關(guān)。

3.2 河冰物候特征變化規(guī)律

3.2.1 河冰物候時(shí)間變化規(guī)律

黃河（內(nèi)蒙古段）5個(gè)子段河冰物候節(jié)點(diǎn)均集中在10～128 d（11月中下旬至翌年2月中下旬，圖7）S3～S4段FUS出 現(xiàn) 最 早，平 均 為14.7 d（11月24日），S1段出現(xiàn)最晚，平均約為38.5 d（12月18日）。2015—2020年各段FUS均有不同程度的偏晚現(xiàn)象，S1和S2段表現(xiàn)較為明顯，分別偏晚約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波動(dòng)范圍較大，波動(dòng)幅度約15 d，其余各段波動(dòng)幅度均在10 d左右。各段FUE主要集中在12月中旬至1月上旬，F(xiàn)UE最早和最晚分別是S3和S1段，近5年平均值分別約為24.6 d（11月24日）和48.2 d（12月28日）；除S2段，F(xiàn)UE波動(dòng)范圍均超過20 d，其余各段均在10 d左右。各段從FUS至FUE所用時(shí)間不盡相同，S1～S5段近5年平均約為9.7、14.0、8.2、12.8 d及10.8 d，其中S2段2016—2017年用時(shí)最長(zhǎng)，為36 d。

圖7 試驗(yàn)段關(guān)鍵物候節(jié)點(diǎn)時(shí)間分布Fig.7 Time distribution diagram of key phenological nodes in the test section

相比較而言，S1～S5段BUS較為分散，可從83 d（S1段）開始一直持續(xù)到126 d（S4段），大多子段的BUS相繼出現(xiàn)在110～125 d，各段之間差異性較大，S1與S4段之間較為明顯，兩段BUS相差約22 d；除S4段，其余各段BUS均有不同偏早現(xiàn)象，變化速率分別為0.8、2.6、2.6 d·a-1及7.0 d·a-1。S1段近5年BUS波動(dòng)范圍較大，約27 d；S4段波動(dòng)范圍穩(wěn)定在5 d以內(nèi)，其余三段波動(dòng)范圍在15 d左右。各段BUE主要集中在三月，其中S1段最早（113 d），S4段最晚（129 d）。與BUS變化規(guī)律相似，除S1和S4段，其余三段BUE均有不同程度偏早現(xiàn)象，變化速率分別為1.6、1.0 d·a-1及7.2 d·a-1。在研究時(shí)段內(nèi)，5段由BUE至FUS的平均封凍期為77 d，其中S1段封凍期最短（47 d），S4段封凍期最長(zhǎng)（96 d）。

進(jìn)一步分析S1～S5段物候信息發(fā)現(xiàn)，S1和S4段從BUS至BUE的 橫向凍 結(jié)速 率 正以0.2 d·a-1和1.4 d·a-1的速度逐年增加，而S3和S5段的橫向凍結(jié)速率呈1.4 d·a-1和1.0 d·a-1的速度逐年減少，S2段近5年基本保持不變，由于每年FUS均集中發(fā)生于S3～S5段，表明黃河內(nèi)蒙古段近5年首封位置逐步向下游偏移。S1～S3段的橫向消融速率呈增加趨勢(shì)，分別為0.4、1.0、1.8 d·a-1，S4段基本保持不變，而S5段以0.6 d·a-1橫向消融速度呈放緩趨勢(shì)。

3.2.2 河冰物候凍結(jié)（消融）空間模式

圖8～9給出了2015—2020年試驗(yàn)段物候期FUS及BUE的凍結(jié)（消融）空間模式?？傮w來講，研究時(shí)間內(nèi)5段間FUS和BUE存在較大差異；對(duì)于FUS，一般S1和S2段最先流凌，然后逐步向其他各段延伸，總體呈“中間早，兩邊晚”的發(fā)展模式；FUE的發(fā)展模式與FUS相反，大體呈“兩邊早，中間晚”的格局，表明黃河（內(nèi)蒙古段）封開河逆序且各段逐步發(fā)展的空間模式。以S1～S5段進(jìn)行分析，S1段河冰FUS發(fā)展模式為“兩岸淺水區(qū)—庫區(qū)中心—庫尾深水區(qū)”，而BUE則呈“庫區(qū)中心—兩岸淺水區(qū)及深水區(qū)”的模式發(fā)展；S2～S5段相似，發(fā)展模式為“兩岸或河灘—主流區(qū)”；需要說明的是，彎道處FUS的發(fā)展模式為“凹岸—凸岸—主流區(qū)”，而BUE則按照“主流區(qū)—凸岸—凹岸”的順序逐步發(fā)展。

圖8 試驗(yàn)段FUS時(shí)空分布Fig.8 Temporal and spatial distribution of FUS in test section

對(duì)S1～S5段進(jìn)一步分析發(fā)現(xiàn)，F(xiàn)US和BUE的空間分布可以反映出水深及河底地形地貌的差異性。以S1段為例，張磊等［37］的研究表明，該河段的河底地勢(shì)由東南向西北傾斜，深水區(qū)位于距壩體約860 m處的庫尾段，淺水區(qū)位于兩岸及西南部，結(jié)合S1段FUS可知，河冰的FUS空間模式與河底地形地貌分布具有一致性，而BUE的空間分布與FUS恰好相反。

4 討論

冬季，由于河道上空霧、云及霾的存在，光學(xué)數(shù)據(jù)在識(shí)別河冰物候信息中受到限制，合成孔徑雷達(dá)數(shù)據(jù)在物候識(shí)別方面效果較好［33］。同時(shí)，黃河河冰物候研究受限于數(shù)據(jù)的時(shí)相分辨率，因此，常在研究中出現(xiàn)某一物候期缺失數(shù)據(jù)的情況［18］。與大尺度海冰（湖冰）物候監(jiān)測(cè)和植被物候監(jiān)測(cè)不同，本文基于Sentinel-1時(shí)間序列數(shù)據(jù)，綜合曲線斜率法和動(dòng)態(tài)閾值法構(gòu)建了黃河（內(nèi)蒙古段）5個(gè)子段的河冰物候曲線，提出了S1～S5段河冰關(guān)鍵物候節(jié)點(diǎn)的最優(yōu)提取指標(biāo)，解決了河冰物候信息難以提取的問題［38］。本文發(fā)現(xiàn)，時(shí)序數(shù)據(jù)提取、擬合以及提取方法是影響提取精度的重要因素。在時(shí)間序列后向散射系數(shù)提取中，由于河道清溝和濕雪的存在，易將河冰誤提取成水體，可能導(dǎo)致物候提取的BIAS較大，故應(yīng)選用穩(wěn)封期影像進(jìn)行控制點(diǎn)的選取［25，34］。在數(shù)據(jù)擬合時(shí)，與S-G濾波、Lowess平滑等重構(gòu)河冰物候曲線的方法相比，Logistic擬合更能反映河冰生消規(guī)律，同時(shí)不受時(shí)序數(shù)據(jù)等間距的限制，適用于黃河冰情解譯［14，24］。各河段數(shù)據(jù)量不同，但擬合效果良好，主要與樣本類型及特征有關(guān)［36，39］。在提取方法方面，由于各段后向散射系數(shù)的差異，單一閾值易造成FUS和BUE提取的失真，動(dòng)態(tài)閾值根據(jù)各段物候曲線特性提取物候期，但無法量化導(dǎo)致各段之間差異性的不確定性因素（如含沙量、反射率、太陽高度角、河底地質(zhì)構(gòu)造等）。

黃河（內(nèi)蒙古段）各子河段物候特征與前人湖冰物候特征的研究相比，既具統(tǒng)一性，也具差異性。S1段是典型庫區(qū)型河段，水流流速相對(duì)較小，F(xiàn)US和BUE的發(fā)展模式與湖冰發(fā)展模式基本相同［20］。S2～S5段順直部分的FUS和BUE發(fā)展模式與S1段相似，區(qū)別在于，彎道處由于水流離心力的存在導(dǎo)致河冰生長(zhǎng)多集中于凹岸側(cè)，其FUS發(fā)展模式為“凹岸—凸岸—主流區(qū)”。同時(shí)，彎道凹岸處由于存在大量堆冰，使得河岸灘在凍融侵蝕和淘刷作用下易崩塌，故彎道凹岸側(cè)應(yīng)加強(qiáng)防護(hù)［40-42］。

黃河河冰物候特征不僅受氣溫、風(fēng)速、太陽輻射以及積雪覆蓋等的影響，并且也受水位、流量，河道邊界條件以及上下游水庫防凌調(diào)度的影響［43-45］，其中氣溫和水力條件被認(rèn)為是影響黃河河冰物候信息的主要因子［46］。例如：氣溫越高初冰日日期越晚；水位越高，河道與大氣表面接觸面積就越大，河冰凍結(jié)所需要的熱量條件也隨之增大，從而導(dǎo)致河流消融時(shí)間變得緩慢等［47］。因此，在未來研究中，如何在模型中顯示結(jié)合氣象和水文數(shù)據(jù)是一個(gè)值得考慮的方向，這可能會(huì)進(jìn)一步提高河冰表面物候的提取精度。

圖9 試驗(yàn)段BUE時(shí)空分布Fig.9 Temporal and spatial distribution of BUE in test section

5 結(jié)論

本文基于2015—2020年間Sentinel-1/2雷達(dá)/光學(xué)影像數(shù)據(jù)，利用Logistic函數(shù)擬合，綜合曲線斜率法及動(dòng)態(tài)閾值法確定了5個(gè)試驗(yàn)段河冰關(guān)鍵物候節(jié)點(diǎn)的閾值，分析了河冰物候特征變化規(guī)律及其主要驅(qū)動(dòng)因素，得到如下結(jié)論：

（1）基于Sentinel-1時(shí)間序列數(shù)據(jù)，綜合斜率法和動(dòng)態(tài)閾值法可用于黃河（內(nèi)蒙古段）河冰物候提取，整體提取精度較高。

（2）S1～S5段FUS和BUE的河冰物候最優(yōu)閾值提取指標(biāo)分別為：0.1 Sa、0.2 Sa、0.1 Sa、0.05 Sa以及0.05 Sa，BIAS在±3 d以內(nèi)；斜率法提取FUE和BUS的BIAS在±5 d以內(nèi)。

（3）S1～S5段近5年FUS平均變化速率分別為1.4、1.0、0.8、0.2 d·a-1及0.4 d·a-1；S2段 近5年FUS波動(dòng)范圍較大，波動(dòng)幅度為15 d左右，其余各段波動(dòng)幅度均在10 d左右；除S4段，其余各段BUS均有不同程度偏早現(xiàn)象，變化速率分別為0.8、2.6、2.6 d·a-1及7 d·a-1；S1段近5年BUS波動(dòng)范圍較大，約為27 d；S4段波動(dòng)范圍穩(wěn)定在5 d以內(nèi)，其余三段波動(dòng)范圍在15 d左右。

（4）S1段FUS發(fā)展模式為“兩岸淺水區(qū)—庫區(qū)中心—庫尾深水區(qū)”；S2～S5段FUS發(fā)展模式為“兩岸或河灘—主流區(qū)”；彎道處FUS發(fā)展模式為“凹岸—凸岸—主流區(qū)”；S1段BUE發(fā)展模式為“庫區(qū)中心—兩岸淺水區(qū)及深水區(qū)”；S2～S5段順直河段BUE發(fā)展模式為“主流區(qū)—兩岸或河灘”；彎道處BUE發(fā)展模式為“主流區(qū)—凸岸—凹岸”。