基于星载雷达特征参数的黄河冰厚反演及冰储量估算

刘 斌，冀鸿兰，翟涌光，张宝森，郜国明

（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）×104m³增加到（2 437±412）×104m³；河冰完全冻结期，冰储量无较大变化；河冰消融期，气温开始回 暖，2020年2月8日至3月15日冰储量从（2 238±413）×104m³减 小 到（474±289）×104m³。2016—2021年5个凌汛期的冰储量变化如图11所示，各年凌汛期内冰储量随时间呈先增加，后保持相对稳定，最后减小的整体变化特征，这与凌汛期的气温变化是相符的，2016—2021年各年度冰储量逐年增加，各年度稳封期平均冰储量分别为1 070×104、1 882×104、2 124×104、2 335×104、2 500×104m³，结合图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）×104m³增加到（2 437±412）×104m³；河冰完全冻结期，冰储量无较大变化；河冰消融期，气温开始回暖，2020年2月8日至3月15日冰储量从（2 238±413）×104m³减小到（474±289）×104m³。

（5）2016—2021年各年度内冰储量随时间呈先增加，后保持相对稳定，最后减小的整体变化特征，2016—2021年的年际冰储量逐年增加。