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

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

0 引言

淡水冰（例如冰川、湖冰、河冰）是冰凍圈的重要組成部分，在全球氣候、地表面能量平衡、水分通量、水濱及水生生態(tài)系統(tǒng)健康中扮演重要角色。在人口聚居區(qū)域，其對(duì)社會(huì)、經(jīng)濟(jì)層面的影響同樣不容小覷。河水凍結(jié)導(dǎo)致的冰凌災(zāi)害作為中高緯度河流易生自然災(zāi)害，嚴(yán)重威脅沿岸人民生命財(cái)產(chǎn)安全，制約沿岸可持續(xù)發(fā)展。以黃河內(nèi)蒙古段為例，每年冰期可持續(xù)4～5個(gè)月，受河道形態(tài)、水文條件、氣象條件、人類活動(dòng)影響，每年凌情不盡相同［1］，河段內(nèi)冰凌災(zāi)害頻發(fā)［2］。據(jù)已公開文獻(xiàn)記載［3-6］，1993—2008年間所發(fā)生的冰凌災(zāi)害中累計(jì)受災(zāi)群眾達(dá)2.4×104人，受災(zāi)面積接近240 km2，保守估算經(jīng)濟(jì)損失超11×108元人民幣。然而，相比冰凍圈其他組成部分，對(duì)于河冰的認(rèn)知仍存在差距［7］，河冰空間分布變化事關(guān)堤防的鞏固安全、上游水庫防凌調(diào)度以及跨河建筑物的布設(shè)等防凌減災(zāi)工作，因此了解河冰的空間分布是十分必要的。

科研人員［8-10］對(duì)河冰變化研究多集中于基于水文觀測(cè)的封開河時(shí)間、封凍時(shí)長(zhǎng)、槽蓄水增量以及凌峰流量等凌情數(shù)據(jù)的變化。但是，隨著氣候、人類活動(dòng)、河道特征、水沙特性等因素變化影響，河冰的空間分布同樣會(huì)發(fā)生變化，而空間分布的變化勢(shì)必會(huì)影響防凌減災(zāi)工作的對(duì)策。黃河內(nèi)蒙古段封凍河段總長(zhǎng)度大于700 km，若進(jìn)行實(shí)地河冰分布觀測(cè)需耗費(fèi)大量人力、物力、財(cái)力、時(shí)間，且無法保證觀測(cè)標(biāo)準(zhǔn)的統(tǒng)一性和數(shù)據(jù)的時(shí)效性。1980年以來，遙感的蓬勃發(fā)展為淡水冰研究提供了可靠的技術(shù)手段和大量的歷史遙感影像。目前，基于光學(xué)遙感技術(shù)的淡水冰研究總體呈現(xiàn)三個(gè)方向：其一，季節(jié)性凍結(jié)冰的物候研究。例如，Kropá?ek等［11］利用中分辨率成像光譜儀（MODIS）8天合成數(shù)據(jù)反演青藏高原59個(gè)湖泊2001—2010年間湖冰的凍融日期；姚曉軍等［12］利用MODIS及Landsat TM/ETM+遙感影像反演并分析可可西里地區(qū)湖冰物候變化，得到了湖冰凍融的空間模式為一岸擴(kuò)展向另一岸；Chu等［13］利用紅外波段在冰和水表面的不同反射率，對(duì)2000—2015年間加拿大西北領(lǐng)地奴河的封開河時(shí)間進(jìn)行了估計(jì)，精度令人滿意。其二，多年凍結(jié)冰的空間變化研究。例如，在冰川監(jiān)測(cè)領(lǐng)域上，劉娟等［14］利用Landsat OLI遙感影像結(jié)合冰川編目數(shù)據(jù)，得出了岡底斯山冰川面積萎縮加快及北朝向萎縮最快的結(jié)論；李志杰等［15］利用Landsat衛(wèi)星影像，對(duì)1993—2016年間喀喇昆侖山什約克流域冰川空間變化進(jìn)行分析，得出冰川整體萎縮且北朝向萎縮速率最快的結(jié)論；高永鵬等［16］利用Landsat OLI影像對(duì)河西內(nèi)流區(qū)冰川變化進(jìn)行監(jiān)測(cè)。其三，冰凌監(jiān)測(cè)研究。K??b等［17-18］利用高時(shí)間、空間分辨率遙感影像分別對(duì)西伯利亞地區(qū)Lena River、Amur River及北美Yellowstone River、Yukon River浮冰流速及冰塞進(jìn)行監(jiān)測(cè)，所得結(jié)果具有較高精度。趙水霞等［19］使用Landsat8衛(wèi)星數(shù)據(jù)對(duì)黃河什四份子段河冰生消及冰塞過程進(jìn)行了分析；李超等［20］利用Landsat8衛(wèi)星7、4、3波段組合獲得良好的河冰提取精度，為冰情解譯提供最佳波段組合；楊中華等［21］提出“四星三源”（MODIS、ETM+、CBERS-02、RADARSAT）遙感動(dòng)態(tài)凌汛監(jiān)測(cè)模式，實(shí)現(xiàn)黃河凌情實(shí)時(shí)監(jiān)測(cè)。

綜上，受上游水庫調(diào)度、河道萎縮等因素影響，作為季節(jié)性凍結(jié)的河冰同樣存在時(shí)空變化。因此，為補(bǔ)充黃河內(nèi)蒙古段河冰分布資料，判定河冰偎堤高危區(qū)，預(yù)警凌汛險(xiǎn)情，本文基于歷史遙感影像，使用遙感、GIS技術(shù)提取河冰分布信息并分析河冰時(shí)空分布及變化特征，以期為黃河內(nèi)蒙古段重點(diǎn)堤段巡查、防凌調(diào)度、跨河建筑物布置等提供科學(xué)依據(jù)。

1 研究區(qū)概況

黃河內(nèi)蒙古段起點(diǎn)為內(nèi)蒙古自治區(qū)烏海市，流經(jīng)烏蘭布和沙漠東緣，穿過陰山山脈南麓與鄂爾多斯高原之間的沖積平原，終到呂梁山西側(cè)進(jìn)入山西境內(nèi)（圖1）。流域內(nèi)屬溫帶大陸性氣候，全年少雨，冬春處于冷高壓控制之下，每年11月末至12月初寒潮或冷空氣入侵致使干流開始流凌、封凍，翌年3月氣溫回暖河道開始解凍、河冰消融，冰期歷時(shí)3～4個(gè)月。封河時(shí)，河冰由下游向上游逐漸封凍，開河時(shí)順序相反由上游開始解凍。研究河段沿程經(jīng)過五市一盟共17個(gè)旗縣。截至2017年末，沿岸17個(gè)旗縣常住人口452.16×104人，占全自治區(qū)人口17.71%；糧食總產(chǎn)量487.66×104t，占全自治區(qū)糧食總產(chǎn)量的15%［22］；京呼銀蘭通訊光纜、河套灌區(qū)、西北電網(wǎng)高壓輸電線路、達(dá)拉特電廠、托克托電廠、包鋼分布在河段兩岸［23］。研究河段全長(zhǎng)約644 km，起點(diǎn)為海勃灣水庫壩址終至托克托縣河口村（頭道拐水文站附近）。為便于分析，按水文站位置、河道比降、河型，將研究區(qū)分為4個(gè)子段［24］，各子段特征見表1。

圖1 研究區(qū)地理位置及遙感影像覆蓋范圍Fig.1 Geographic location of the study area and the boundary of remote sensing image

表1 研究區(qū)分段及各子段特征Table 1 sub-reaches of the study area and characteristics of different sub-reaches

2 數(shù)據(jù)來源及研究方法

2.1 數(shù)據(jù)選取及預(yù)處理

本研究收集1989—2019年間覆蓋研究區(qū)的所有影像，成像時(shí)間盡可能保證為每年1—2月，篩選出云量小于10%且河冰不被云覆蓋的影像數(shù)據(jù)，如表2所示。其中91景Landsat5、7、8衛(wèi)星所攜帶的TM、ETM+、OLI傳感器影像（后文簡(jiǎn)稱L5TM、L7ETM+、L8OLI）及2景Sentinel-2衛(wèi)星MSI傳感器（后文簡(jiǎn)稱S2MSI）影像，因2019年研究區(qū)中R3、R4段Landsat8衛(wèi)星影像質(zhì)量不高，故選用相近日期的Sentinel-2衛(wèi)星影像替代。其中，Landsat系列衛(wèi)星數(shù)據(jù)獲取自美國(guó)地質(zhì)調(diào)查局（United States Geological Survey，https：//earthexplorer.usgs.gov/），Landsat

表2 遙感影像數(shù)據(jù)Table 2 Remote sensing data

衛(wèi)星數(shù)據(jù)為經(jīng)過地面控制點(diǎn)校正的Level 1Tie級(jí)數(shù)據(jù)，覆蓋研究區(qū)的影像行號(hào)Row為32，條帶號(hào)Path分別為127、128、129；Sentinel-2衛(wèi)星數(shù)據(jù)獲取自歐洲航天局（European Space Agency，https：//scihub.copernicus.eu/dhus/#/home），Sentinel-2衛(wèi)星數(shù)據(jù)為L(zhǎng)2A級(jí)大氣底層反射率影像，覆蓋R3、R4段影像圖幅號(hào)為T49TDE、T49TCE。Landsat系列衛(wèi)星數(shù)據(jù)為經(jīng)過地面控制點(diǎn)幾何校正的影像數(shù)據(jù)，其預(yù)處理過程為：首先，基于遙感圖像處理軟件ENVI中的輻射校正工具對(duì)傳感器進(jìn)行輻射定標(biāo)，將圖像DN值轉(zhuǎn)換為反射率值；其次，利用FLAASH工具進(jìn)行大氣校正處理，獲得地物表面反射率；最后，裁剪圖像減少非研究區(qū)地物對(duì)分類過程的影響，提高數(shù)據(jù)處理效率。其中2012年、2013年所使用L7ETM+數(shù)據(jù)經(jīng)過條帶修復(fù)，所提取河冰信息僅供參考。Sentinel-2衛(wèi)星的L2A級(jí)數(shù)據(jù)在發(fā)布前已經(jīng)過輻射校正、幾何校正等預(yù)處理，可直接使用。

2.2 分類方法

2.2.1 河冰信息提取方法

目前基于光學(xué)遙感數(shù)據(jù)提取冰雪信息的主要方法為指數(shù)法，原理為根據(jù)冰雪與其他地物光譜特征的差異設(shè)置閾值提取冰雪信息。目前歸一化雪指數(shù)（normalized difference snow index，NDSI）被廣泛應(yīng)用在湖冰、積雪、冰川等研究領(lǐng)域，是目前最常用且精度較高的方法［25］。NDSI由Dozier于1989年首次提出［26］，Riggs等［27］基于NDSI提出SOMAP算法，Hall等［28］精確了SNOMAP算法提取雪信息的閾值。SNOMAP算法提取冰雪信息需滿足如下條件：

式中：NDSI=(ρGreen-ρSWIR1)·(ρGreen+ρSWIR1)-1，ρGreen為綠光波段反射率值；ρSWIR1為短波紅外波段反射率值；ρNIR為近紅外波段反射率值。需要指出的是NDSI閾 值 取 值 范 圍 在0.25～0.45之 間，常 用 值為0.4。

通過使用NDSI閾值提取的冰雪信息中包含一定量的水體信息，SNOMAP算法中設(shè)置近紅外波段反射率閾值的目的就是剝離水體信息。鑒于河冰中存在清溝（未封凍水體），本文使用并改進(jìn)SNOMAP算法進(jìn)行河冰信息提取。因黃河河冰內(nèi)含沙量較大、冰內(nèi)氣泡較多，加之河冰表面粗糙程度不同，勢(shì)必會(huì)影響河冰光譜特征［29］，若利用SNOMAP算法首先需對(duì)NDSI提取冰水混合類的閾值進(jìn)行試驗(yàn)。由于Hall對(duì)NDSI剝離除冰、水的其他地物像元的效果已經(jīng)做了詳盡的描述，在此不再贅述。試驗(yàn)主要以土壤、河冰、水體作為主要地物類型，根據(jù)各種地物的NDSI指數(shù)的分布情況確定合理閾值。參照參考影像分別選取三種地物樣本，導(dǎo)出樣本點(diǎn)像元在各波段反射率值并計(jì)算NDSI，衛(wèi)星遙感影像與參考影像信息如表3所示。如圖2所示，土壤像元的NDSI值多為負(fù)值，僅有ETM+影像中3個(gè)像元的NDSI指數(shù)值大于0；河冰與水體的NDSI值大體保持在0.25以上，僅有4個(gè)河冰像元NDSI值（2005年2個(gè)，2016年1個(gè)，2019年1個(gè)）在0.2左右。結(jié)合NDSI的常用取值范圍并在盡量排除異常值的情況下，選取0.25作為提取冰水混合類的閾值取值。此外，對(duì)比河冰與水體像元NDSI指數(shù)值的箱型圖可知，二者NDSI指數(shù)值有重疊部分，說明單純使用NDSI閾值無法準(zhǔn)確提取河冰像元，需其他方法區(qū)分兩種地物像元。

圖2 主要地物純像元NDSI指數(shù)值分布圖Fig.2 NDSI value distributions of pure pixels of major land cover types

表3 衛(wèi)星遙感影像及參考影像介紹Table 3 Introduction of satellite remote sensing image and reference image

2.2.2 水體信息提取方法

由于NDSI與NIR波段反射率值的閾值組合提取河冰像元出現(xiàn)了水體像元被誤分為河冰的現(xiàn)象，故提出歸一化未封凍水體指數(shù)（NDUWI）與NDSI的閾值組合提取水體信息方法以獲得較高精度的河冰像元。因河道內(nèi)未封凍水體對(duì)藍(lán)、綠光反射率差距較大，而河冰對(duì)藍(lán)、綠光反射率差距較小，因此，結(jié)合河冰及水體的光譜特征差異，提出歸一化未封凍水體指數(shù)（Normalized Difference Unfrozen Water Index，簡(jiǎn) 稱NDUWI）。NDUWI計(jì) 算 公 式如下：

式中：ρBlue為藍(lán)光波段反射率值。

為驗(yàn)證NDUWI提取效果，根據(jù)河冰、水體像元的近紅外波段（NIR）反射率值、歸一化水體指數(shù)（NDWI）值、改進(jìn)的歸一化水體指數(shù)（MNDWI）值及歸一化未封凍水體指數(shù)（NDUWI）值，確定其分布情況，再根據(jù)分布情況確定區(qū)分二者的閾值。NDWI、MNDWI計(jì)算公式如下：

式中：ρNIR為近紅外波段反射率值。

如圖3（a）所示，NIR波段反射率值在TM及MSI傳感器影像中提取效果較好，但TM傳感器在NIR波段區(qū)分兩種地物的閾值在0.2左右，而非SNOMAP算法中0.11，ETM+、OLI影像中，河冰、水體NIR反射率值并無明顯的閾值；如圖3（b）、3（c）所示，河冰、水體的NDWI和MNDWI（NDSI）指數(shù)值之間無明顯的閾值存在（其中MNDWI計(jì)算公式與NDSI一致，故圖中MNDWI指數(shù)值分布圖相當(dāng)于放大NDSI指數(shù)值分布圖）；如圖3（d）所示，河冰、水體像元的NDUWI值之間的存在清晰的閾值，ETM+、MSI影像中閾值為0.2，TM、OLI影像中閾值為0.25。NIR反射率值及3種指數(shù)提取水體像元閾值見表4。

表4 NIR反射率值、NDWI、MNDWI、NDUW閾值Table 4 Threshold for NIR reflectance，NDWI，MNDWI and NDUWI

圖3 河冰、水體純像元NIR波段反射率值及不同指數(shù)值分布圖Fig.3 Different indexes and NIR reflectance value distributions of pure pixels of river ice and water cover types

以2016年2月1日OLI影像為例展示水體提取效果，當(dāng)ρNIR大于0.11時(shí)，透明度較高的河冰像元被誤分為水體，隨著閾值的降低仍然出現(xiàn)誤分現(xiàn)象，當(dāng)閾值低至0.045以下時(shí)，同時(shí)出現(xiàn)河冰像元誤分、水體像元漏分的現(xiàn)象［圖4（b）］。NDWI、MNDWI提取水體信息結(jié)果相似，均出現(xiàn)河冰被誤分水體現(xiàn)象［圖4（c）、4（d）］。如圖4（e）所示，利用NDUWI閾值提取水體像元效果較為理想。4個(gè)傳感器影像利用NDSI與NDUWI閾值組合提取河冰、水體的效果如圖5所示。

圖4 4種方法提取水體效果對(duì)比Fig.4 Contrast of extraction result by 4 methods

圖5 利用NDSI與NDUWI閾值組合提取河冰、水體信息效果Fig.5 Extraction results of 4 scenes from 4 sensors classified by NDSI and NDUWI

2.2.3 精度驗(yàn)證

用于分類驗(yàn)證的真實(shí)地物輪廓是基于Google Earth高清歷史影像人工目視解譯的。由于高清歷史影像的空間分辨率在0.454 m，目視解譯的誤差忽略不計(jì)。選取2017年1月9日巴彥淖爾市附近Google Earth歷史影像進(jìn)行人工目視解譯，驗(yàn)證同日成像的Landsat影像的分類精度。驗(yàn)證過程將使用高清影像進(jìn)行人工解譯的河冰、水、土壤作為真值，計(jì)算混淆矩陣。然后分別計(jì)算各類生產(chǎn)者、用戶者、總體精度及Kappa系數(shù)。如表5所示，土壤分類精度最高，冰與水的分類精度相近?？傮w精度可達(dá)89.55%，Kappa系數(shù)為0.88，分類結(jié)果較為理想。

表5 冰、水、土壤分類精度評(píng)估Table 5 Evaluation of classification accuracy in river ice，water and soil

2.3 空間分析方法

得到河冰分類信息后生產(chǎn)河冰矢量數(shù)據(jù)，并對(duì)矢量數(shù)據(jù)進(jìn)行分布變化分析。在進(jìn)行分布變化分析時(shí)，R1段選取右岸國(guó)家干線公路（G110線）作為參照物（因R1段兩岸堤防級(jí)別較低，遙感影像中不易辨認(rèn)），R2～R4段以兩岸堤防為參照物。R1段將G110線作為起點(diǎn)建立每級(jí)寬度為0.5 km的12級(jí)緩沖區(qū)，R2、R3段以兩岸堤防作為起點(diǎn)分別建立每級(jí)寬度為0.5 km的5級(jí)緩沖區(qū)，R4段以R2、R3的方法建立4級(jí)緩沖區(qū)。計(jì)算河冰面矢量與緩沖區(qū)重疊部分的面積，分析每級(jí)緩沖區(qū)內(nèi)河冰面積變化，進(jìn)而得到河冰分布變化。河冰面積計(jì)算方法假設(shè)構(gòu)成矢量多邊形的點(diǎn)坐標(biāo)為（xi，yi），則多邊形的面積S為：

其中，按順時(shí)針依次選點(diǎn)計(jì)算多邊形面積為正值，逆時(shí)針依次選點(diǎn)計(jì)算多邊形面積為負(fù)值。

3 研究結(jié)果

3.1 黃河內(nèi)蒙古段河冰總體分布

1989—2019年間黃河內(nèi)蒙古段穩(wěn)封期河冰主要分布在中游（圖6），中游兩岸堤防之間大部分河灘被河冰覆蓋，特別是在1999—2015年間，中游兩岸堤防之間少有裸露河灘，有河冰偎堤現(xiàn)象。上游段2015年前河冰分布較為穩(wěn)定，無顯著變化，2015年起上游首部河冰消失，其余河段河冰寬度收窄。下游段河冰寬度較小，漫灘現(xiàn)象較輕，多數(shù)年份可分辨主槽位置。

圖6 1989—2019年黃河內(nèi)蒙古段河冰分布Fig.6 Spatial distribution of river ice of Yellow River in Inner Mongolia section from 1989 to 2019

3.2 河冰面積變化

3.2.1 黃河內(nèi)蒙古段河冰面積變化

如圖7所示，1989—2019年黃河上游內(nèi)蒙古段河冰大致經(jīng)歷“穩(wěn)定期”（1989—1997）、“擴(kuò)張期”（1998—2000年）、“萎縮期”（2001—2019年）。2000年河冰擴(kuò)張至最大范圍，全段河冰面積為1 723.50 km2；2019年萎縮至最小，面積為936.27 km2；多年平均值為1 349.54 km2（由于冰期跨越相鄰兩年，后文均使用翌年年份代表冰期年份，例如始于1999年終于2000年的冰期以2000年表示）。

1989—1997年間，河冰面積在1 072.93 km2至1 380.32 km2之間波動(dòng)，平均值為1 223.68 km2，河冰面積維持在較低水平；1998—2000年間，河冰經(jīng)歷快速上升階段，3年間累積增加面積為561.91 km2，所擴(kuò)張面積占2000年（31年間最大值）河冰面積的32.6%；2001—2015年間，河冰緩慢萎縮，年平均萎縮率為-18.4 km2·a-1，其中2003年出現(xiàn)河冰面積極小值，其原因?yàn)榉夂悠冢?002年12月上旬）黃河內(nèi)蒙古段上游的石嘴山站流量較往年少23.2%且氣溫為往年最低，封河期氣溫低、流量小、水位低，使河道內(nèi)快速形成穩(wěn)定的冰蓋。穩(wěn)封后流量保持在200 m3·s-1，使河水在冰蓋下，故河冰面積出現(xiàn)極小值。2016年河冰面積驟減472.07 km2，同比減少32.64%；2017年、2018年出現(xiàn)小幅回升，2019年出現(xiàn)最小值，3年間河冰面積在947.38 km2至1 177.56 km2之間波動(dòng)，平均值為1 048.47 km2，低于31年平均值，河冰面積處于低水平。

3.2.2 各子河段河冰面積變化

如圖8所示，除1992年、2003年出現(xiàn)極小值外，對(duì)應(yīng)表1中的R1段1989—2014年河冰面積穩(wěn)定，多年 均 值為189.27 km2；2015年、2016年 分別萎 縮66.16 km2、66.93 km2，2016年河冰面積僅為2014年40.88%。R2段大體呈現(xiàn)“先增后減”的變化特征，1989—2008年間河冰面積波動(dòng)上升，其中1997—1999年連續(xù)增加177.82 km2，2000—2008年河冰面積保持高位并于2008年達(dá)到31年最大值（801.49 km2）；2009年起波動(dòng)下降，其中2016年由771.75 km2驟降至581.12 km2，2017年、2018年小幅回升，2019年回落至552.04 km2，與1989年水平相近。R3、R4段河冰面積變化特征大體與R2段一致，呈現(xiàn)為“先增后減”，但波動(dòng)幅度較R2段劇烈。1996—1999年經(jīng)歷河冰面積大幅增加時(shí)期，其中僅1999年兩段同比增加了51.69%、72.84%，河冰面積達(dá)到歷史最大值；2006—2019年河冰面積波動(dòng)下落，僅2016年兩段同比下降46.72%、32.60%，2017年兩段表現(xiàn)略有差異，2017年、2018年R3段出現(xiàn)小幅回升、R4段則保持穩(wěn)定。但與R2段不同的是，R3、R4段在2001—2005年間出現(xiàn)明顯“下降—回升”的過程。

圖8 1989—2019黃河內(nèi)蒙古段各子河段河冰面積變化曲線Fig.8 Variations of river ice area in each sub-reach in Inner Mongolia section of Yellow River from 1989 to 2019

3.3 河冰空間分布變化

R1段分布變化主要表現(xiàn)為沿流向的縱向萎縮。由圖9（b）、9（c）可知，R1段左右岸除1992年、2003年有較大波動(dòng)外，2015年前各緩沖區(qū)河冰面積較為穩(wěn)定；2015年、2016年各緩沖區(qū)河冰面積同時(shí)大幅減少［見圖9（a）中2015年河冰分布］，2018年稍有回升。對(duì)比兩岸邊緣區(qū)，經(jīng)過2015年、2016年的兩次河冰萎縮過程，左岸4.5～5.5 km兩個(gè)緩沖區(qū)河冰基本消失，右岸0～0.5 km緩沖區(qū)河冰已經(jīng)消失，由此可見，發(fā)生縱向萎縮同時(shí)，R1段也出現(xiàn)河冰橫向萎縮，且左岸橫向萎縮更為顯著。

圖9 R1段河冰分布圖及兩岸各緩沖區(qū)河冰面積年際變化Fig.9 Distribution of river ice in section R1 and interannual variation of river ice area in each section along both banks

R2段河冰變化以沿流向橫向分布變化為主。如圖10（a）所示，1989年河冰主要以河道中心分布，而2019年河冰以右岸分布為主。對(duì)比圖10（b）、10（c）可知，左右岸分布變化主要出現(xiàn)在靠近堤防的邊緣區(qū)（0～1 km），1989—1998年左右兩岸邊緣河冰面積呈現(xiàn)波動(dòng)上升趨勢(shì)，左岸外緣波動(dòng)幅度較右岸波動(dòng)幅度大；1999年兩岸河冰擴(kuò)張，左右岸邊緣河冰面積分別同比增加64.84%、50.84%；2000—2010年河冰面積保持小范圍波動(dòng)，在此期間，兩岸河冰大致呈現(xiàn)對(duì)稱分布；2011年起左右岸同時(shí)開始萎縮，左岸邊緣區(qū)河冰年均萎縮率為2.17 km2a-1、2.48 km2·a-1，右岸為1.14 km2·a-1、0.88 km2·a-1，且右岸0～1 km兩個(gè)緩沖區(qū)面積始終大于左岸相應(yīng)區(qū)域。至2019年左岸0～1 km緩沖區(qū)河冰面積基本恢復(fù)1989年水平，其余緩沖區(qū)低于1989年水平；而右岸各緩沖區(qū)保持或高于1989年水平。由圖10（a）可知，R2段河冰漫灘現(xiàn)象嚴(yán)重，2000年兩岸堤防間基本被河冰覆蓋，裸露河灘面積較小，且河冰偎堤現(xiàn)象嚴(yán)重。隨時(shí)間推移，2016年、2019年左岸河冰萎縮，出現(xiàn)大面積裸露河灘；而右岸仍有河冰堆積，且河冰距離堤防距離較近，部分位置仍然存在偎堤現(xiàn)象。

圖10 R2段河冰分布圖及兩岸各緩沖區(qū)河冰面積年際變化Fig.10 Distribution of river ice in section R2 and interannual variation of river ice area in each section along both banks

R3段分布變化集中在左右岸邊緣區(qū)，31年間經(jīng)歷向兩岸擴(kuò)張?jiān)傧蛑鞑畚s的過程［圖11（a）］，2016年以后大體以主槽河冰為主。如圖11（b）、11（c）所示，1989—2000年間，左岸邊緣區(qū)河冰面積小于左岸其他緩沖區(qū)河冰面積；2000—2010年間除2003年、2004年外，邊緣區(qū)外側(cè)（0～0.5 km）河冰面積已超過其他緩沖區(qū)內(nèi)河冰面積，但與其他各緩沖區(qū)河冰面積相差不大；2000年后各緩沖區(qū)基本保持同步變化。右岸各緩沖區(qū)變化趨勢(shì)大體一致，邊緣區(qū)河冰面積變化幅度最大且始終小于其他各緩沖區(qū)，各緩沖區(qū)河冰面積隨距堤距離的增加而增加。對(duì)比圖11（a）各代表年河冰分布位置可知，1989年河冰以主槽為軸線向兩岸擴(kuò)散，2000—2015年間河冰漫灘現(xiàn)象嚴(yán)重，以致左右岸河冰均有偎堤現(xiàn)象出現(xiàn)。2016年起河冰萎縮，漫灘偎堤現(xiàn)象大為緩解；2019年基本以主槽封凍為主。

如圖12（b）、12（c）所示，R4段左右岸各緩沖區(qū)河冰面積年際變化較大，2019年除0～0.5 km緩沖區(qū)回歸1989年水平外，其他各緩沖區(qū)河冰面積均萎縮。由圖12（a）可知，R4段首部河冰漫灘偎堤現(xiàn)象最為嚴(yán)重，而河段尾部河冰分布較為穩(wěn)定，其余部分以局部漫灘為主。相比2015年前各年份，2016年起漫灘現(xiàn)象有較大緩解，裸露河灘面積增大；2019年河冰主要以主槽封凍為主。

4 討論

4.1 海勃灣水庫對(duì)下游冰情影響

2014年前海勃灣壩址下游河道常年封凍，2014年5月海勃灣水庫運(yùn)行后，R1段由于封凍長(zhǎng)度的縮短導(dǎo)致其河冰面積萎縮幅度最大，冬季壩址至巴彥木仁段不再封凍。原因是：水庫運(yùn)行后，庫內(nèi)蓄水溫度呈現(xiàn)分層現(xiàn)象，庫底水溫保持恒定，導(dǎo)致下泄水流“冬暖夏涼”［30］進(jìn)而影響下游河道水溫；其次，水庫的泄洪閘改變了水流的邊界條件，使水流內(nèi)部急劇摻混，流速的改變導(dǎo)致下泄至河道內(nèi)的水體在一定距離內(nèi)不封凍。對(duì)于巴彥高勒站以下河段（R2、R3段），封河期海勃灣水庫控制下泄流量在650 m3·s-1，對(duì)入庫流量實(shí)行“多蓄少補(bǔ)”［31］，水庫對(duì)流量的平穩(wěn)控制有效避免流量波動(dòng)，使冰花順利流向下游，有效緩解下游冰情。以2014年末河冰封凍過程為例，海勃灣水庫投入運(yùn)行后，理論上河冰面積會(huì)有所下降，然而下游各河段均出現(xiàn)河冰面積的大幅上升，造成該現(xiàn)象的原因?yàn)?014年冬季封河期流量為847 m3·s-1，超過預(yù)案設(shè)定下泄流量197 m3·s-1，比多年平均值多300 m3·s-1；2014年11月25日海勃灣水庫下泄流量由0 m3·s-1激增至1 420 m3·s-1，至11月27日下泄流量仍保持在960 m3·s-1［32］，下泄水流流入下游河段時(shí)遭遇寒流入侵，致使河水快速凍結(jié)后壅冰上岸，造成河冰面積增加。由圖7可知，下游R2、R3、R4段受大流量過程影響，2015年水庫投入運(yùn)營(yíng)后河冰面積并未下降反而增加，而2016年起河冰面積開始下降并在之后處于低水平；從圖11（a）可知，2015年R4段上、中河冰漫灘現(xiàn)象較為嚴(yán)重，2016年后得到明顯緩解?？梢姡２獮乘畮炜捎行Ь徑庀掠瘟枨?，且水庫平穩(wěn)調(diào)度對(duì)下游防汛工作至關(guān)重要，科學(xué)的防凌調(diào)度可有效控制河冰擴(kuò)散，緩解下游河冰漫灘現(xiàn)象。

4.2 河床演變對(duì)河冰分布變化的影響

河冰過程包括了復(fù)雜的相互作用，包括水力學(xué)、力學(xué)以及熱力學(xué)，同時(shí)也受天氣和水文條件的影響［33］，而河床正是水力學(xué)當(dāng)中的邊界條件，河床作為河水的承載體，而河水又是河冰的承載體，其性質(zhì)勢(shì)必會(huì)影響河冰分布情況。主槽位置變化、深度變化、寬度變化、灘槽高差、彎道數(shù)量等特征都影響著河道內(nèi)河水結(jié)冰封凍的過程。黃河上游內(nèi)蒙古段存在兩個(gè)典型河段，游蕩型的巴彥高勒至三湖河口段以及彎曲型的包頭至頭道拐段。R2段河冰面積為四段中最大，而R4段河冰分布變化屬四段最顯著。

R2段淤積最為嚴(yán)重、演變最為復(fù)雜［34］。河床寬度、主槽位置不穩(wěn)定，為河冰分布造成諸多不確定性。河槽淤積與河床變窄共同作用造成灘槽落差減小，容易造成河水漫灘、輸冰不暢的局面。安催花等［35］研究表明，截至2012年R2段、R3段、R4段河槽淤積嚴(yán)重，斷面萎縮；梅艷國(guó)等［36］研究表明，R2段在1977—2000年河床寬度總體呈減小趨勢(shì)，其中1985—1995年間河床寬度呈現(xiàn)波動(dòng)下降；2000—2014年保持在小范圍波動(dòng)但總體不變。這意味著河道向窄而淺的斷面形態(tài)發(fā)展。相應(yīng)地，由圖8中R2段河冰面積變化曲線可知，R2段河冰面積從1989年的593.32 km2增加至2000年的779.39 km2，增幅為16.92 km2·a-1，與河床寬度減小趨勢(shì)相對(duì)應(yīng)，其中1989—1995年河冰面積呈現(xiàn)波動(dòng)上升與河床寬度波動(dòng)下降相對(duì)應(yīng)；2000—2014年河冰面積保持656.66～801.49 km2之間波動(dòng)（其中僅有3年低于720 km2分別為2003年、2011年、2014年），與河床寬度保持穩(wěn)定相對(duì)應(yīng)。由此可見，河道淤積后河底高程抬升、河床寬度變窄導(dǎo)致河床橫斷面向窄而淺發(fā)展，過水面積的減小容易導(dǎo)致冰凌堵塞，從而造成河冰堆積壅水漫灘進(jìn)而造成河冰面積的增加。

卡冰現(xiàn)象的發(fā)生是由河道地形、流量、冰產(chǎn)量等因素共同作用產(chǎn)生。泥沙淤積造成河床的抬高導(dǎo)致主槽冬季輸冰時(shí)濕周減小、過流面積減小，一旦上游來冰量劇烈增加，極易在河道收窄、彎道處形成冰塞。上游海勃灣水庫若按照防凌調(diào)度計(jì)劃保證流量平穩(wěn)，流凌期、封凍期大流量穩(wěn)封期小流量，可在主槽中形成穩(wěn)定冰蓋，減小河冰卡冰漫灘（河冰漫灘減少的表現(xiàn)形式就是河冰面積的減少）。穩(wěn)封后控制下泄流量可使水流在冰蓋與主槽形成的類似管道的邊界條件下輸水輸冰，不會(huì)造成水鼓或冰上過流從而造成河冰的擴(kuò)散。

河道水流的平面形式和發(fā)展趨勢(shì)稱為河勢(shì)，其中涵蓋河道水流動(dòng)力軸線位置、走向以及河灣、岸線和沙洲、心灘等分布與變化的趨勢(shì)。R4段屬?gòu)澢秃佣?，幾乎不存在沙洲、心灘，河?shì)變化主要以水流走向、河彎、岸線變化為主。如圖13所示，對(duì)比主槽變化圖發(fā)現(xiàn)，R4段主槽擺動(dòng)主要發(fā)生在公布營(yíng)子至白家圪旦、九股地至東李三壕內(nèi)，其余河段河道主槽變化并不明顯。對(duì)比主槽變化，公布營(yíng)子至侉圪堵之間變動(dòng)較大，相比于1988年流凌前，1990年、1999年、2005年、2008年、2016年均發(fā)生主槽擺動(dòng)，截至2019年公布營(yíng)子至侉圪堵段增加5處彎道，侉圪堵彎道附近河道逐漸由左岸擺動(dòng)至右岸；九股地至東李三壕之間1990年、1995年、1999年、2003年、2005年、2016年、2018年均發(fā)生擺動(dòng)彎道，九股地至東李三壕段增加4個(gè)彎道，東李三壕附近彎道主要以萎縮為主。對(duì)比1989年以來河冰分布情況（圖12），公布營(yíng)子、侉圪堵、九股地以及東李三壕幾處河道擺動(dòng)頻繁的區(qū)域常有堆冰現(xiàn)象發(fā)生。東李三壕至頭道拐水文站之間河段主槽31年間擺動(dòng)幅度非常小，河冰堆積情況相比其他河段較輕?？梢?，河勢(shì)變化較大河段河冰漫灘現(xiàn)象比較嚴(yán)重。

圖12 R4段河冰分布圖及兩岸各緩沖區(qū)河冰面積年際變化Fig.12 Distribution of river ice in section R4 and interannual variation of river ice area in each section along both banks

圖13 1989—2019年R4段主槽位置變化Fig.13 Changes of main channel location in R4 reach from 1989 to 2019

5 結(jié)論

基 于1989—2019年 間Landsat衛(wèi) 星TM、ETM+、OLI傳感器及Sentinel-2衛(wèi)星MSI傳感器影像數(shù)據(jù)，提出了歸一化未封凍水體指數(shù)（NDUWI），利用NDUWI指數(shù)值與NDSI指數(shù)值的閾值相結(jié)合的方法，提取了黃河河冰信息并進(jìn)行時(shí)空變化分析，主要結(jié)論如下：

（1）黃河上游內(nèi)蒙古段河冰分布呈現(xiàn)“中間多兩邊少”，巴彥高勒至三湖河口河冰面積最大，包頭至頭道拐最小。研究區(qū)內(nèi)河冰經(jīng)歷“平穩(wěn)期”（1989—1997年）、“擴(kuò)張期”（1998—2000年）、“萎縮期”（2001—2019年），全 段 河 冰 面 積 在936.27～1723.50 km2之間，河冰面積最大、最小值分別出現(xiàn)在2000年和2019年。

（2）1989—2014年海勃灣水庫壩址至巴彥高勒之間（R1段）河冰分布穩(wěn)定，受海勃灣水庫建成使用的影響，2015年起河冰面積開始下降并保持低水平。巴彥高勒至三湖河口（R2）、三湖河口至包頭（R3）、包頭至頭道拐（R4）河冰面積變化趨勢(shì)大體保持一致，三段均于1998—2000年呈現(xiàn)上升趨勢(shì)，2001年起呈現(xiàn)下降趨勢(shì)，其中2016年出現(xiàn)明顯下降后并保持低水平。

（3）各子河段中，R1段河冰以沿流向萎縮為主，2015年后海勃灣水庫壩址下游20 km內(nèi)不再封凍。R2段2001年前左右兩岸大體呈現(xiàn)對(duì)稱變化，2001年起左岸邊緣區(qū)（0～0.5 km、0.5～1 km）年平均萎縮率為2.17 km2·a-1、2.48 km2·a-1，右岸為1.14 km2·a-1、0.88 km2·a-1，左岸河冰萎縮而右岸大體保持不變。R3段河冰分布呈現(xiàn)向兩岸擴(kuò)張?jiān)傥s的變化特征，2016年前漫灘現(xiàn)象嚴(yán)重，2016年起兩岸同時(shí)向主槽萎縮，漫灘現(xiàn)象大為緩解。R4段2016年前河段首部漫灘現(xiàn)象嚴(yán)重，尾部以主槽封凍為主，其余河段局部漫灘，2016年起以全段主槽封凍為主。

（4）海勃灣水庫合理控制下泄流量可有效地減少黃河內(nèi)蒙古段河冰漫灘現(xiàn)象；建議著重加強(qiáng)巴彥高勒至三湖河口段右岸堤防巡查。