基于星載雷達特征參數的黃河冰厚反演及冰儲量估算

劉 斌，冀鴻蘭，翟涌光，張寶森，郜國明

（1.內蒙古農業大學水利與土木建筑工程學院，內蒙古 呼和浩特 010018；2.黃河水利委員會黃河水利科學研究院，河南 鄭州 450003）

0 引言

冰凌洪水是高緯度河流所特有的一種自然災害，嚴重威脅著河道水工建筑物的安全與河岸生態系統的穩定［1-2］。冰厚是冰凌成因分析、冰情模擬及預報的重要基礎信息，可為防凌減災提供重要依據［3-4］；無論是冰的防治，還是冰的利用，冰厚均是一個關鍵的參數，也是比較難以監測的物理指標，如何對其進行準確、有效地估算一直是河冰研究中的重點和難點［5］。黃河內蒙古段每年凌汛期長達4個月，河道蜿蜒曲折，是凌災重點防控河段，其冰厚信息的獲取對黃河防凌減災具有重要意義［6］。

目前，現場人工鉆孔測量冰厚是最為可靠的測量手段，但受自然環境和交通條件等因素的限制，耗時費力，成本較高，且具有一定的危險性［7］。度日法冰厚估算模型自Stefan［8］提出以來，在冰工程領域得到廣泛應用［9-10］，國內學者王軍等［11］對該模型進行改進，用以估算黃河內蒙古段的河冰冰層厚度，該方法雖能較好地反映河冰的生消趨勢，卻難以體現冰厚的空間分布。近年來，探地雷達［12-14］、機載雷達［15］在黃河冰厚監測中的應用，極大地促進了河冰研究的發展，但仍難以進行大范圍的河冰厚度時空監測。遙感技術的飛速發展，使得以非接觸手段進行大范圍冰厚監測成為可能［16］，遙感估算冰厚主要采用光學遙感、被動微波遙感和主動微波遙感相結合等方法開展［17］。光學遙感法一般基于海冰厚度反照率反演模型，在渤海海冰研究中取得了一系列進展［18-21］，由于河冰的類型比海冰復雜得多，不同厚度的河冰光譜曲線規律性較差，且光學遙感數據易受天氣影響，限制了光學遙感法在河冰厚度反演中的應用。被動微波遙感法一般基于不同頻率的亮溫差來進行冰厚的反演，由于當前被動微波遙感數據空間分辨率較低，多應用在海冰或湖冰厚度的研究中［22-24］。主動微波遙感法一般通過分析星載雷達參數與冰厚的相關性，進而建立經驗模型反演冰厚，由于主動微波遙感數據擁有媲美光學遙感數據的空間分辨率，且不受天氣影響，有利于進行全天候、全天時的冰情監測，在加拿大河流的冰厚反演中應用較多［25-27］，但國內相關研究極少，因此探究主動微波遙感法在黃河冰厚反演中的可行性尤為迫切。

Sentinel-1遙感數據作為目前使用最廣泛的免費星載合成孔徑雷達（SAR）數據，對河冰厚度的反演有巨大潛力。本文基于黃河凌汛期野外冰厚實測數據，分析Sentinel-1雷達特征參數與河冰厚度之間的相關性，并建立冰厚反演經驗模型，最后對黃河凌汛期冰厚進行反演，分析冰厚的時空分布特征，并對冰儲量進行估算，以期為黃河內蒙古段防凌減災提供參考依據。

1 研究區及數據來源

1.1 研究區概括

黃河內蒙古段地處黃河流域最北端，介于39°～41°N，106°～112°E之間，海拔在1 000 m左右，2016—2021年5個凌汛期平均氣溫為-6.5℃，依次流經石嘴山、巴彥高勒、三湖河口、包頭及頭道拐5個水文站。因其地理位置特殊、氣象條件復雜，導致凌汛期封開河逆序。凌汛期一般從11月下旬開始，到翌年3月中下旬結束，歷時100余天。本研究區范圍選取包頭與頭道拐水文站之間（圖1），該河段為彎曲型河段，易發生卡冰結壩現象［28］。

圖1 研究區概況Fig.1 Overview map of the study area

1.2 數據來源

1.2.1 野外實測數據

野外實測地點位于頭道拐水文站附近的什四份子彎道處，該河段穩封期河冰冰層較厚，便于野外踏勘，在此進行連續三年的河冰厚度信息采集，采集日期分別為2019年1月9日至13日、2020年1月16日至20日、2021年1月12日至16日。依照中華人民共和國水利行業標準《河流冰情觀測規范》進行斷面布置，勘測過程中每隔10～20 m布設1個采樣點，進行人工鉆孔測量冰厚，共選取15個斷面合計300多個采樣點（圖1中斷面按年份順序分別用紅、藍、綠三種顏色進行標注），出于安全考慮，未對清溝附近冰厚進行人工測量。利用RTK采集冰厚采樣點的地理坐標，便于后期與遙感影像進行空間匹配；利用量冰尺測量采樣點的河冰厚度，有積雪的同時測量雪深。

1.2.2 衛星遙感數據

Sentinel-1是歐空局進行的哥白尼計劃中的第一個任務，該任務是由A、B兩顆衛星在同一軌道平面組成的星座。A星和B星分別于2014年4月和2016年4月發射，單星重返周期為12天，雙星組合重返周期為6天，均載有波長約5.6 cm的C波段合成孔徑雷達，可提供全天時、全天候的高分辨率雙極化產品。Sentinel-1雷達影像包含強度與極化信息，用于分析星載雷達特征參數與河冰厚度的相關性。

Sentinel-2是歐空局進行的哥白尼計劃中的第二個任務，該任務目前也是由A、B兩顆衛星在同一軌道平面組成的星座。A星和B星分別于2015年6月和2017年3月發射，單星重返周期為10天，雙星組合重返周期為5天，均攜帶多光譜成像儀，可覆蓋13個光譜波段，影像幅寬為290 km，空間分辨率最高可達10 m。各波段詳細信息見表1。Sentinel-2光學影像在本文用于提取黃河的主河道。

表1 Sentinel-2各波段詳細參數Table 1 Sentinel-2 detailed parameters of each waveband

本文選取研究區2016—2021年5個凌汛期的57景Sentinel-1雷達影像及對應精密軌道數據、5景Sentinel-2光 學 影 像。5景Sentinel-2光 學 影 像 均 在凌汛期前采集，采集時間分別為2016年10月27日、2017年10月31日、2018年10月26日、2019年11月10日、2020年11月4日。衛星遙感數據均來源于歐空局官網（European Space Agency，ESA，https：//scihub.copernicus.eu/dhus/#/home/）。

1.2.3 氣象數據

研究區的凌汛期日均氣溫數據來自黃河網（http：//www.yrcc.gov.cn/），降雨和降雪數據來自托克托縣氣象站。氣象數據的詳細信息見圖2，以2019—2020年凌汛期為例。已有研究表明，雷達特征參數對氣溫、降水變化較為敏感，故需獲取氣溫和降水數據以保證研究的準確性［27］。

圖2 2019—2020年黃河內蒙古段凌汛期氣象信息及雷達影像獲取時間Fig.2 Meteorological information and radar image acquisition time during the ice flood season of the Inner Mongolia section of the Yellow River from 2019 to 2020

河冰生消過程與氣溫息息相關，尤其是累積凍冰度日（cumulative freezing degree day，CFDD）與累積融冰度日（cumulative thawing degree day，CTDD）等參數，CFDD和CTDD分別為日均氣溫低于和高于冰點的度數在一定時間范圍的累積［29］，計算公式定義如下：

式中：CFDD和CTDD分別是累積凍冰度日與累積融冰度日（℃·d）；θ為累積凍冰度日與3倍累積融冰度日的差值（℃·d）；Ta為研究區日均氣溫（℃）；Tf為冰凍結和融化的臨界溫度（℃），本文采用-5℃［30-31］；ds和de分別為黃河內蒙古段凌汛期的開始和結束日期。

2 研究方法

2.1 基于Sentinel-2的河道邊界提取

2.1.1 預處理

大氣校正的目的是消除大氣和光照等因素對地物反射的影響，獲得地物真實的表面反射率。從歐空局獲取的Sentinel-2影像L1C級數據是已完成輻射校正與幾何校正的大氣頂層（TOA，top of atmosphere）反射率產品，需對其進行大氣校正得到大氣底層（BOA，bottom of atmosphere）反射率產品，以表征地物真實的光譜特性。SNAP（Sentinel application platform）軟件的大氣校正插件sen2cor結合了用于執行大氣校正的先進技術，這些技術是針對Sentinel-2制定的，可對L1C級數據進行大氣、地形和卷云校正，得到L2A級數據。

2.1.2 黃河主河道提取

為獲取黃河主河道邊界，需剔除陸地區域，消除陸地對河冰信息提取的干擾。利用改進的歸一化差異水體指數（MNDWI）可以突出光學影像中的水體，且能夠抑制植被信息［32］，計算公式如下［33］：

式中：Green為綠色波段的像素值；SWIR為短波紅外波段的像素值。

然后通過目視解譯選擇合適的閾值，使水體與陸地區分明顯，對經過處理的影像進行二值化及矢量化處理。將初步獲取的河道矢量與Sentinel-2影像中的真實水陸分界線對比，對分界線進行精細化修改，誤差控制在一個像素以內，即小于10 m，最終獲得黃河主河道邊界。

2.2 基于Sentinel-1的河冰厚度反演及驗證

2.2.1 預處理

（1）后向散射系數提取

采用SNAP軟件的流程圖工具對Sentinel-1影像的SLC數據依次進行軌道校正、去熱噪聲、輻射定標、Deburst、多視、極化濾波、地形校正及分貝化處理，最終獲得VV與VH兩種極化方式的后向散射系數。地形校正采用美國宇航局（National Aeronautics and Space Administration，NASA，http：//srtm.csi.cgiar.org/srtmdata/）發布的空間分辨率為30 m的SRTM DEM數據。

（2）極化特征參數提取

利用PolSARpro軟件對Sentinel-1影像的SLC數據依次進行軌道校正、輻射定標、Deburst、生成極化矩陣、地形校正、極化濾波及極化分解處理，獲取4種極化特征參數：α、極化熵、λ1、λ2。α描述了探測目標內部的自由度，說明了目標平均散射的“類型”［34］。極化熵較好地反映了目標散射的隨機性［35］。λi表示極化相干矩陣的第i個特征值［36］。

2.2.2 冰厚反演

本研究采用統計回歸方法建立冰厚反演模型，主要包括3個步驟，分別是：相關性分析、模型建立、模型驗證。首先繪制雷達特征參數與實測冰厚的散點圖，并計算雷達特征參數與實測冰厚之間的Pearson相關系數，分析雷達特征參數與冰厚的相關性，公式如下：

式中：r為Pearson相關系數；xˉ為變量x的均值；yˉ為變量y的均值。

選擇與冰厚相關性最高的雷達特征參數建立線性回歸模型，從經時空匹配的野外實測冰厚數據中隨機選取2/3作為模型構建數據，用于確定模型的待定系數，其余1/3作為模型的驗證數據，采用均方根誤差RMSE、平均相對誤差MRE、平均絕對誤差MAE等表征模型精度。

2.2.3 冰儲量估算

根據反演的河冰厚度可對研究區的冰儲量進行估算，如公式（6）：

式中：V為冰儲量（m3）；A為遙感影像中單個像元表示的地面面積（m2）；hi為遙感影像所反演的冰厚（m）。

3 結果與分析

3.1 黃河主河道提取結果

研究區2016—2021年5個凌汛期前的黃河主河道如圖3所示，該河段為彎曲型河段，彎多且彎曲系數大，主河道邊界年際變化顯著。整體來看，凌汛期前主河道寬度逐年增大，2016—2017年變化最為顯著；彎道數量呈減小趨勢，河道擺動幅度減小并逐步趨于穩定。

圖3 2016—2021年5個凌汛期前的黃河主河道Fig.3 The main channels of the Yellow River before the 5 ice flood seasons from 2016 to 2021

3.2 冰厚反演模型建立及驗證

以往的研究通常采用歸一化技術來限制衛星軌道不同所引起的雷達特征參數差異［37］，本文在每個凌汛期均采用同一軌道的影像進行研究。雷達特征參數與實測冰厚的散點圖如圖4所示，2020年各雷達特征參數與實測冰厚的一致性均較差，分布比較散亂。這是由于2020年在冰期野外踏勘期間，試驗區河冰冰層均被較厚的積雪覆蓋，且期間有降雨現象，冰層上方形成了濕雪，影響了各雷達特征參數與冰厚之間的相關性。2019年、2021年在冰期野外踏勘期間，冰層上無積雪。雷達特征參數與實測冰厚的Pearson相關系數見表2，對于無濕雪覆蓋的冰而言，VV極化后向散射系數與河冰厚度的相關性最高，λ1次之。

表2 極化特征參數與實測冰厚的Pearson相關系數Table 2 Pearson correlation coefficients of polarization characteristic parameters and measured ice thickness

圖4 極化特征參數與實測冰厚的散點圖Fig.4 Scatter plot of polarization characteristic parameters and measured ice thickness

兩種基本的數值模擬方法通常被用來估算冰層厚度：基于冰記錄的隨機技術和基于物理冰生長和消融原理的參數化方法［38］。這兩種模型都需要大量的參數輸入，主要是環境數據［39］，卻都無法提供冰厚的空間分布信息。本文選擇與河冰厚度相關性最高的VV極化后向散射系數建立冰厚反演模型，在雷達特征參數與河冰厚度的相關性分析中，已將SAR影像的各特征參數與實測冰厚數據進行了時空匹配，基于Sentinel-1 VV極化后向散射系數和野外實測冰厚的反演模型如圖5所示，得到公式（7）：

圖5 冰厚反演模型Fig.5 Ice thickness inversion model

冰層的后向散射系數每增加1 dB，冰層的厚度增加5.036 cm，所擬合線性模型的調整R2值為0.657，適用于10～100 cm的冰厚反演，10 cm以下的冰厚認為是水體或薄冰。圖6表明驗證數據的均方根誤差RMSE為9.82 cm，平均相對誤差MRE為13.46%，平均絕對誤差MAE為8.26 cm。

圖6 冰厚反演結果Fig.6 Ice thickness inversion results

3.3 冰厚反演結果

圖7顯示了2019—2020年黃河凌汛期包頭至頭道拐水文站河段的冰厚反演結果，河冰厚度平均值見于表3。研究區的河冰平均厚度范圍為6.9～58.2 cm，2019年11月16日為凌汛期前3天，此時河道以水體為主，彎道處存在少量薄冰；2019年11月中下旬至2020年1月初為河冰生長期，河冰平均厚度從6.9 cm增加到58.2 cm，清溝逐步減小；2020年1月初至2月初為河冰完全凍結期，河冰平均厚度基本穩定在55 cm左右；2020年2月初至3月中旬為河冰消融期，河冰平均厚度從53.2 cm減小到16.1 cm。表3顯示，在河冰生長與消融期，氣溫參數θ與河冰厚度的增減一致，較好地反映了河冰的生消變化，但是氣溫參數θ無法判斷河冰完全凍結期。

表3 2019—2020年黃河內蒙古段凌汛期冰厚與氣溫參數θTable 3 Ice thickness and air temperature parameters θ of the Inner Mongolia section of the Yellow River during the ice flood season from 2019 to 2020

圖7 2019—2020年黃河凌汛期冰厚反演結果Fig.7 Ice thickness inversion results of the Yellow River during the ice flood season from 2019 to 2020

圖8顯示了2019—2020年黃河凌汛期什四份子彎道處的冰厚反演結果。河冰生長期，彎道處首先發生凍結，凹岸早于凸岸；河冰完全凍結期，彎道處冰厚分布基本穩定，凸岸冰厚相對較小，冰厚在30～50 cm居多，凹岸冰厚相對較大，冰厚在70～80 cm居多，清溝大小略有變化；河冰消融期，凸岸融化早于凹岸。冰層下測速后發現彎道處的主流在凸岸側，流速較大，因此冰底沖刷侵蝕嚴重，水流帶來的浮冰和冰花無法平穩堆積；由于彎道獨特的地形，大量冰盤、浮冰及冰花向凹岸移動并堆積，且凹岸流速較小冰盤并置，造成凹岸冰層較厚。

圖8 2019—2020年黃河凌汛期什四份子彎道冰厚反演結果Fig.8 Ice thickness inversion results of Shisifenzi Bend of the Yellow River during the ice flood season from 2019 to 2020

圖9顯示了2016—2021年5個黃河穩封期的平均冰厚、平均封凍面積及平均氣溫，2016—2019年各穩封期平均冰厚呈增大趨勢，分別為40.4 cm、49.3 cm、49.7 cm、55.7 cm，同時期平均封凍面積也持續增大，應該是由于同時期平均氣溫逐年下降所致；2019—2021年穩封期由于平均氣溫上升，導致平均冰厚減小，而平均封凍面積不降反增，可能是由于該河段河勢、水力等條件變化所致。

3.4 冰儲量估算結果

本文冰儲量通過河冰厚度與封凍面積的乘積計算，冰儲量估算的誤差可通過誤差傳播公式確定，河冰厚度反演的均方根誤差Δh為9.82 cm，封凍面積A的誤差主要由冰水分類的精度確定，此處不做考慮，因此冰儲量估算的中誤差為ΔV=A·Δh=9.82×10-2Am3，下面對年內冰儲量與年際冰儲量進行分析。

圖10顯示了2019—2020年黃河凌汛期包頭至頭道拐水文站河段的冰儲量，河冰生長期，隨氣溫降低，2019年11月16日至2020年1月3日冰儲量從（280±162）×104m3增加到（2 437±412）×104m3；河冰完全凍結期，冰儲量無較大變化；河冰消融期，氣溫開始回 暖，2020年2月8日至3月15日冰儲量從（2 238±413）×104m3減 小 到（474±289）×104m3。2016—2021年5個凌汛期的冰儲量變化如圖11所示，各年凌汛期內冰儲量隨時間呈先增加，后保持相對穩定，最后減小的整體變化特征，這與凌汛期的氣溫變化是相符的，2016—2021年各年度冰儲量逐年增加，各年度穩封期平均冰儲量分別為1 070×104、1 882×104、2 124×104、2 335×104、2 500×104m3，結合圖9可以發現2016—2020年各穩封期平均氣溫持續下降，與同時期平均冰儲量變化相符；2019—2021年穩封期平均氣溫升高、平均冰厚減小，而平均封凍面積增大，導致同時期平均冰儲量不僅沒有減小，反而有所增加。

圖9 2016—2021年5個黃河穩封期平均冰厚、平均封凍面積及平均氣溫Fig.9 Average ice thickness，average frozen area and average temperature of the Yellow River during the five frozen periods from 2016 to 2021

圖10 2019—2020年凌汛期冰儲量Fig.10 Ice storage during the ice flood season from 2019 to 2020

圖11 2016—2021年5個凌汛期的冰儲量Fig.11 Ice storage of 5 ice flood seasons from 2016 to 2021

4 討論

一般來說，雷達后向散射系數的大小取決于成像目標的介電特性、散射特性和傳感器特性［40］。

材料的介電特性是由其介電常數表征的。對于淡水冰來說，這個常數幾乎不隨冰溫而變化［41］，與海冰不同［42］，但覆蓋雪的含水量對其影響較大［43］。表4總結了淡水冰研究中幾種重要物質的介電常數［40，44-47］。電磁波在介質中的傳播過程所引起的電磁損耗由介質的介電常數決定。當介質的介電常數較大時，電磁損耗高，穿透深度小，例如，液態水是一種高損耗介質，當雪或冰中的水分含量很高時，電磁波的穿透深度就會減小，對于5 GHz的C波段電磁波而言，積雪中僅5%的體積含水量就會使穿透深度減小到10 cm左右［48］。相反，低介電常數意味著低水分含量，例如，淡水冰是一種低損耗介質，頻率在1～10 GHz下的電磁波，對淡水冰的穿透深度從100 m到10 m［48］，這給微波探測冰物候及冰厚變化提供了可能性［49］。

表4 淡水冰研究中重要物質的介電常數Table 4 Dielectric constants of important substances in freshwater ice research

河冰的散射特性主要由兩方面決定：表面散射和體積散射。表面散射以鏡面反射和漫反射為主。鏡面反射發生在較為光滑的冰表面，在微波的入射波長尺度內，可認為是光滑的。光滑潮濕的冰面（即冰層上有積水）和平靜的水面，都會遠離微波入射方向處鏡面反射入射微波的能量，因此后向散射響應較弱，雷達圖像中顯示較暗。相反，漫反射發生在相對粗糙的表面，在所有方向上幾乎均勻的反射微波，并將其中一部分入射能量反射回傳感器，因此一般粗糙的冰面在雷達圖像中會顯得相對明亮。當入射波穿透冰面或干燥的薄雪，并在冰層內的介電不連續處重復反射時，就會發生體積散射，這些不連續處物質的尺寸需與入射微波的波長相似［43］，河冰中影響體散射的主要因素有：裂縫、氣泡和雜質，張邀丹等［50］對黃河冰內的氣泡已有較為深入的研究。也有研究表明，微波在管狀氣泡的前向散射與冰水界面的高介電對比度相結合，會發生一種“雙反射機制”，進而產生較強的雷達回波［51］，而河冰的電磁輻射傳輸建模研究發現冰水界面對雷達后向回波貢獻較大［52］。

合成孔徑雷達的傳感器特性主要由其波長、極化方式和入射角決定。前兩者對于傳感器來說是固定的，這里主要討論入射角的影響。星載雷達通過改變不同的微波入射角進而以短于重返周期的時間間隔對同一地理區域成像，增強其動態監測的潛力。已有研究發現，同一目標的后向散射信號隨其微波入射角的變化而變化［53］。

本文僅從星載雷達特征參數與河冰厚度的相關性出發，探究河冰厚度反演的有效手段，反演結果雖較好地反映了河冰生消變化，但河冰空間分布的合理與準確性仍有待商榷。關于本次研究，以下幾點需要重點思考：（1）出于野外工作的安全考慮，冬季對河冰厚度信息進行采集時，河冰已生長了一段時間，厚度基本在30 cm以上，而以采樣數據所建立的冰厚反演模型是否也符合河冰生長初期冰厚的變化規律。（2）河冰類型復雜多樣，其內部結構也相差甚遠，采用單一方法的河冰厚度反演可能會在部分類型的河冰上出現較大誤差。如果先對河冰類型進行判斷、再賦以經驗厚度值，結合模型反演結果綜合判斷河冰厚度，有助于進一步提升冰厚反演精度。（3）河冰厚度的變化與流量、水位、氣溫、降水等因素息息相關，將相關水文、氣象因素與遙感模型相結合，有助于更好地了解河冰生消機理。（4）目前，對河冰散射機制的研究還不夠完善，這也限制了冰厚反演精度的提高。河冰生消比海冰、湖冰復雜得多，受水動力條件、熱力學條件、河道形態等多因素共同控制，需進一步開展試驗研究與理論分析，為黃河冰凌探測提供科學的理論依據。

5 結論

本文基于黃河凌汛期野外實測冰厚數據，分析Sentinel-1雷達特征參數與河冰厚度的相關性，并建立冰厚反演模型，估算黃河凌汛期的河冰厚度及冰儲量變化，得到如下結論：

（1）通過分析雷達特征參數與實測冰厚的相關性，發現VV極化后向散射系數與黃河冰厚相關性最高，λ1次之，Pearson相關系數分別為0.81077、0.72409。

（2）經驗證，基于VV極化后向散射系數建立的冰厚反演經驗模型的均方根誤差為9.82 cm，平均相對誤差為13.46%，平均絕對誤差為8.26 cm。

（3）2019—2020年凌汛期研究區河冰平均厚度范圍為6.9～58.2 cm，2019年11月16日為凌汛期前3天，此時河道以水體為主，彎道處存在少量薄冰；2019年11月中下旬至2020年1月初為河冰生長期，河冰平均厚度從6.9 cm增加到58.2 cm，清溝逐步減小；2020年1月初至2月初為河冰完全凍結期，河冰平均厚度基本穩定在55 cm左右；2020年2月初至3月中旬為河冰消融期，河冰平均厚度從53.2 cm減小到16.1 cm。

（4）河冰生長期，隨氣溫降低，2019年11月16日至2020年1月3日冰儲量從（280±162）×104m3增加到（2 437±412）×104m3；河冰完全凍結期，冰儲量無較大變化；河冰消融期，氣溫開始回暖，2020年2月8日至3月15日冰儲量從（2 238±413）×104m3減小到（474±289）×104m3。

（5）2016—2021年各年度內冰儲量隨時間呈先增加，后保持相對穩定，最后減小的整體變化特征，2016—2021年的年際冰儲量逐年增加。