基于MODIS数据的额尔齐斯河上游河冰物候研究

胡婉嫔， 效存德， 谢爱红， 沈永平， 雷华锦

（1.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室，甘肃兰州730000；2.中国科学院大学，北京100049；3.北京师范大学地表过程与资源生态国家重点实验室，北京100875）

0 引言

额尔齐斯河是发源于中国新疆阿尔泰山南坡的一条跨国河流，沿阿尔泰山南麓向西北，经中国、哈萨克斯坦和俄罗斯，在俄罗斯境内流入鄂毕河，后注入北冰洋，是中国唯一流入北冰洋的河流［1-6］。额尔齐斯河干流属于荒原河流，水面宽阔，流速平缓，落差小，能行舟航运［7］。20世纪50至60年代，额尔齐斯河在中苏之间运输建筑材料、石油制品等货物，货运量从630万吨增长至1 850万吨，后因中苏关系恶化，该河航运停止［8］。目前随着全球气候变暖，北极变暖的速度远快于预期，北冰洋冰层融化速度加快，极可能在2050年前后出现夏季无冰状态［9］，北极航道成为各国共同关注的热点，随着中国“一带一路”倡议的稳步推进，“冰上丝绸之路”逐渐走进大众的视野。据报道，设滨海国际运输走廊的邀请，建议双方共同开发和利用海上通道特别是北极航道，打造“冰上丝绸之路”。2018年1月，中国政府发布《中国的北极政策》白皮书，提出中俄共建“冰上丝绸之路”，促进北极地区互联互通和经济社会可持续发展，至此“冰上丝绸之路”得到完整阐释［10-12］。有报道称哈萨克斯坦和俄罗斯政府早在2004年已签订相关协议，清理并拓深哈国边境至鄂木斯克市的河道，使额尔齐斯河实现哈俄之间的跨境通航［13］。若中国与哈萨克斯坦之间的额尔齐斯河上游能够通航，中国货物可以从新疆布尔津县经额尔齐斯河直接运往哈萨克斯坦后，到达北冰洋，运往欧美等国家。据哈萨克斯坦估计，若额尔齐斯河上游通航，不仅可以开发沿途的旅游业，而且水路货运量全年可达120万吨［14］。

北半球许多季节性封冻的河流和湖泊对航运和运输业有很大影响，河湖中浮冰密度太大，冰块太厚或严重封冻，将导致每年内陆水道停运数月［15］。额尔齐斯河为季节性封冻河流，每年10月开始出现岸冰，12月封冻，至次年5月初完全开河，约有120天至160天封冻期。随着“冰上丝绸之路”的建设，额尔齐斯河冰冻期及其变化研究显得尤为必要。

当前，学者们基于遥感数据对河冰及其变化过程进行了大量研究，Chaouch等［16］结合一种基于阈值的决策树图像分类算法，利用MODIS图像探测和绘制萨斯奎哈纳河河冰，并确定河冰范围；Chu等［17］使用MOD09GQ第二波段数据研究加拿大西北地区Slave River河沿岸的河流结冰过程；Unterschultz等［18］利 用RADARSAT-1合 成 孔 径 雷 达（SAR）卫星影像描述河流冰在冬季与崩解期间的特性的可行性。此外，在河湖冰信息提取方面，Sarah等［19］利用MODIS数据研究了河冰开河时间及开河的影响因素，Cai等［20］基于开放水域像素数的阈值提取青藏高原58个湖泊的冻结开始和破裂结束时间。

本文拟采用TIMESAT软件，提取河冰物候信息，并与南湾水文站实测资料、布尔津站与哈巴河站气象观测资料进行对比分析，进而对河冰开封河时间、冰期持续时长及开封河速率等进行分析，以期为关于本河段通航所带来的经济效益研究奠定科学基础，为评估河段航运现实能力及未来潜力提供科学依据。由于资料有限，目前仅研究额尔齐斯河布尔津至斋桑泊河段河冰物候。

1 数据源与研究方法

1.1 研究区概况

额尔齐斯河全长4 248 km，在中国境内546 km，流域面积5.7×104km2，年径流量多达111×108m3。本文研究区为额尔齐斯河流域上游，范围为从布尔津到斋桑泊的额尔齐斯河干流河段（47°40′25.88″～47°48′37.16″N，86°54′40.24″～84°31′12.73″E）（图1），全长约194 km，其中中国境内约115 km，哈萨克斯坦境内约79 km，河流最宽处约350 m，最窄处约100 m。研究区在中国境内属于中国新疆阿勒泰平原地区，地处欧亚大陆中心腹地，纬度较高，属于温带大陆性气候，冬季漫长而寒冷，夏季短促、气温平和，最冷月通常在1月，平均气温为-16℃，最热月在7月，平均气温21℃，地区降水量少，年均降水量131～223 mm，蒸发量大，昼夜温差大，光照充足［21］。

图1 额尔齐斯河研究区Fig.1 Irtysh river research area

1.2 数据集及预处理

MOD09GQ遥感数据源于美国国家航空航天局（NASA）网 站（https：//search.earthdata.nasa.gov/search），是MODIS的反射率产品，时间分辨率为1天，空间分辨率为250 m，提供了两个波段的数据，分别是：第一波段红光波段（620～670 nm）和第二波段近红外波段（841～876 nm）（表1）［22］。研究表明近红外波段对于区别水和雪/冰覆盖最敏感，与红光波段相比，近红外波段更少受到大气条件的影响［23］，故选取第二波段数据作为研究数据，其时间范围为2001—2018年的每年9月1日至5月31日。

表1 MOD09GQ数据说明Table 1 Data specification of MOD09GQ

本文的气象数据选用布尔津站（47°25′12″N，86°31′11″E）及哈巴河站（48°01′48″N，86°14′23″E）2001—2017年的逐日平均气温，用于分析封河期逐日平均气温的变化情况，及其与封河速率的关系。水文数据选用南湾水文站（48°00′N，85°41′E）2001—2018年的河冰开封河时间，用于对比验证遥感数据所得结果。

研究区影像通过谷歌地球（Google Earth）和ArcGIS软件处理得到，使用Google Earth手动矢量化研究区范围（额尔齐斯河布尔津至斋桑泊段）得到研究区矢量图层，然后对MOD09GQ数据进行裁剪，最后提取研究区河流的平均反射率值进行后续研究。

1.3 研究方法

1.3.1 研究路线

本研究将使用HANTS算法去除云噪声，用TIMESAT软件得到河冰物候，具体流程如图3所示。

图3 MODIS的MOD09GQ产品对河冰物候分析总流程Fig.3 General process for river ice phenology analysis of MODIS MOD09GQ

1.3.2 河冰物候定义

河冰物候研究涉及河冰周期性形成与消失，河冰物候受到气候变化以及人类的影响，这与植被物候周期性变化相似，植被物候受季节气候和年际变化的影响，因此植被物候的研究方法亦可用于研究河冰物候。当气温降至0℃时，河水开始冻结，并历经水内冰、薄冰、岸冰、冰覆盖和封冻等阶段［24］。在河流开始结冰之前，MODIS影像的地表反射率有最低值，随着河水开始冻结，反射率逐渐增大。由于光学图像的反射率值取决于地表物体的表面属性而不是内部结构，同样的，河冰反射率仅受河冰表面积雪等的影响，因此，本研究认为当反射率达到最大值时是冰与雪的最大覆盖情况［17］。随着早春气温升高，以及上游开河后流向下游的水流量增加，研究区河段开始开河，直至河流完全无冰，结束开河过程。本文定义通过MODIS探测到河冰出现时为开始封河，当河流无冰时为结束开河，从开始封河到河面最大冰/雪覆盖所经历的时间为封河期，从河面有最大冰/雪覆盖到结束开河所经历的时间为开河期（图2）。

图2 河冰物候定义Fig.2 Phenological definition of river ice

1.3.3 时间序列谐波分析法（HANTS）

MODIS数据具有时间分辨率高、空间覆盖范围广等优点，但其影像常被大量云覆盖，尤其在冬季，部分像元长时间（如10天以上）连续被云遮挡，因此需要采用滤波技术降低噪声，填补MODIS时间序列数据的空白。

时间序列的谐波分析法（HANTS算法）［25-26］是为了去除归一化植被指数（NDVI）时间序列中的云污染值而开发的。HANTS算法的基本思想是计算一个傅立叶级数来对像素级观测的时间序列建模，同时识别相对于时间序列模型的异常值。该算法消除了这些异常值，并用傅里叶级数给出的值将其替换。与快速傅里叶变换（FFT）算法相比，HANTS算法在频率选择和时间序列长度上更具灵活性，此外，由于样本在时间上不要求等距，从时间序列中排除某些点更容易。通过去除时间序列中明显的异常值，提取的谐波比直接的FFT算法可靠得多［27］。HANTS算法已经成功的应用于重建无云的遥感植被数据集（例如归一化植被指数（NDVI）、增强型植被指数（EVI）和叶面积指数（LAI））、重建地表温度数据集以及重建河冰反射率时间序列［27-33，17］。因此，本文采用HANTS算法消除云及其他噪声影响，填补MODIS时间序列数据的空白，重建MODIS时间序列。

1.3.4 TIMESAT

虽然遥感时间序列数据在监测植被季节方面的价值已经得到认可，但从这些序列中探索和提取季节参数的方法有限，因此，TIMESAT软件是为了提取季节参数而开发［34］。现已应用于诸多研究，例如表征物候学［35］、促进生态系统中的数据分类［36］以及监测在火灾季节时人类的足迹［37］等。

该软件实现了三种基于最小二乘拟合的NDVI数据处理方法。第一种方法采用局部多项式函数进行拟合，该方法可分类为自适应Savitzky Golay滤波器（简称S-G滤波）。另外两种方法是普通的最小二乘法，分别为非对称高斯函数和双逻辑函数，其数据适合于不同复杂度的模型函数。这三种处理方法都使用了季节性（单模态或双模态）的初步定义，以及生长季节的近似时间［38-39］。本文使用河冰反射率数据代替NDVI数据，所得生长季节的近似时间即为开封河时间及最大冰/雪覆盖时间。

2 结果与分析

2.1 滤波提取结果对比

S-G滤波器及高斯滤波器提取的河冰物候信息如图4所示，通过比较分析两个时间序列结果，高斯拟合所得时间序列更符合研究所需，最终选用其所得河冰物候信息。由于MODIS部分数据被云遮盖（尤其是冬季影像），所得反射率值高于地物反射率，经过HANTS算法去云以及TIMESAT软件滤波之后，MOD09GQ数据中因为云及其他噪声引起的部分高反射率值得到纠正（图4）。

图4 由MOD09GQ第二波段得到的额尔齐斯河布尔津至斋桑泊段河冰物候Fig.4 River ice phenology of Irtysh river from Burqin to Lake Zaysan section was obtained by MOD09GQ band 2

2.2 河冰物候信息

表2 列出了2001—2018年每年的开始封河时间、结束封河时间、最大冰/雪覆盖时间、最大冰/雪覆盖时所对应的MOD09GQ第二波段的反射率、冰期持续时间、封河速率以及开河速率，其中开封河速率计算方法是开河（或封河）时所对应的河冰反射率与最大冰/雪覆盖时河冰反射率的差值，除以开河期（或封河期）时长。

2.3 结果验证

根据表2得到冰/雪最大覆盖日期，结合预览高分1号影像上研究区的河冰变化情况，选取了2015年的2月2日、3月15日、4月8日以及4月12日的MOD09GQ第二波段的反射率数据进行对比分析。2月2日为2015年冰/雪最大覆盖日期，预览高分1号影像上显示4月8日本文研究区上游有部分冰，而下游斋桑湖附近冰较多，4月12日研究区上游无冰，下游斋桑湖附近无影像。如图5所示，2月2日的反射率值总体最大，之后3月15日、4月8日、4月12日的反射率值明显降低。这表明，随着河冰开河，MOD09GQ的第二波段反射率值逐渐降低，从而认为MOD09GQ的第二波段数据能有效监测河冰变化。

图5 2015年冰/雪最大覆盖期到开河结束期间4天的MOD09GQ第二波段反射率值的变化情况Fig.5 The reflectivity of MOD09GQ second band is worth changing for 4 days from the maximum ice/snow cover period in 2015 to the end of break-up

表2 根据MOD09GQ第2波段数据（841～876 nm）得到的额尔齐斯河布尔津至斋桑泊段2001—2018年间河冰物候信息Table 2 According to MOD09GQ band 2 data（841～876 nm），the river ice phenology information of the Irtysh River from Burqin to Lake Zaysan section from 2001 to 2018 was obtained

由于研究区水文站较少，并且水文站监测河流开封河时间定义与本文不同，因此本文选取额尔齐斯河上，位于布尔津与斋桑泊中间的南湾水文站数据进行对比验证，分别分析由MOD09GQ数据及南湾站得到的2001—2018年间冰期持续时间的变化趋势。由于本文定义MODIS能够探测到河冰出现时为开始冻结，到河流无冰时为开河结束，而水文站定义在封冻期和开河期，当河冰达到一定百分比则认为开始封冻和结束开河，因此本研究所得冰冻期历时大于南湾站数据所得冰冻期历时。2001—2018年间，由MOD09GQ数据及南湾水文站所得冰期持续时间均略有变长趋势［图6（a）、6（b）］，并且通过线性分析得到其相关系数R=0.306，呈正相关［图6（c）］，说明本文所得结果与南湾站实测数据结果的变化趋势一致，均呈现冰期历时延长的趋势。并且水文站及本文研究结果中各年开封河时间变化趋势具有一致性，如2013年南湾站、布尔津站及本文结果解冻日期同比2012年分别提早10天、21天及10天，封冻日期同比2012年分别晚30天、21天及31天，趋势一致，天数相近。

图6 根据MODIS数据所得冰期持续时间与南湾站数据所得冰期持续时间分别的变化趋势及二者的相关关系Fig.6 based on MODIS data from the glacial ice duration and Nanwan station data duration，respectively，the change trend and the correlation of both［Hydrological station in the Nanwan station measured ice duration changes with time（a），by MOD09GQ data of glacial duration changes with time（b），Nanwan station with the duration of the glacial obtained in this paper（c）］

因此，本文认为通过MOD09GQ第二波段数据可以有效监测额尔齐斯河布尔津至斋桑泊段的河冰物候信息。不仅如此，相比于水文站分布稀疏、只能监测较短距离河段物候信息的缺点，遥感数据可以更加灵活地监测不同河段河冰物候信息。然而遥感数据空间分辨率与时间分辨率不能兼得，因此想要监测河冰物候信息，只能选取时间分辨率高而空间分辨率较低的数据，从而遥感数据更适宜于监测大型河流的河冰物候情况。

2.4 开封河时间分析

表2 显示，额尔齐斯河布尔津至斋桑泊段的开始封河时间集中于每年10月，冰/雪最大覆盖时间主要集中于每年1月，结束开河时间主要集中于每年4月。最早开始封河年份为2016年，日期是9月29日，最晚开始封河年份为2013年，日期是11月13日最早封河比最晚封河时间提前45天；最早结束开河在2013年，其日期为3月28日，最晚结束开河为2002年，其日期为4月28日，最早开河比最晚开河时间提前31天，封河与开河时间均提前。分析发现，在2001—2018年这17年中，开始封河时间与结束开河时间总体均有提前趋势［图7（a）、7（b）］，但开始封河时间随年份变化的斜率绝对值（0.623）大于结束开河时间随年份变化斜率的绝对值（0.267）。因此可以认为，封河比开河提前的时间长，导致冰期持续时间将延长［图6（b）］。

图7 开始封河时间与结束开河时间变化趋势（纵轴表示每一年的第几天）Fig.7 The change trend of river freezing time and river breaking time（The vertical axis represents the day of each year），［The change of river break-up time（a），the change of river freeze-up time（b）］

2.5 封河与开河速率变化分析

开河速率与封河速率无明显大小区别，但在速率变化上确有很大不同。通过线性分析，得到封河速率随时间变化的斜率为-0.007，开河速率随时间变化的斜率为0.021，表明2001—2018年，封河速率有变慢趋势，开河速率有增快趋势（如图8）。若额尔齐斯河布尔津段至斋桑泊段通航，则每年秋冬季（即封河期），由于封河速度变慢，河冰形成缓慢，可通航时间将延后。而在每年春季（即开河期），由于开河速度加快，河冰将快速减少，河流达到通航条件时间将提前。因此，随着秋冬季封冻期可通航时间延后，春季开河期可通航时间提前，每年的总通航时间将增加，可带来更大经济效益。

图8 开封河速率变化趋势Fig.8 The change trend of river break-up and river freeze-up rates is that the river freeze-up rate decreases gradually（a）and the river break-up rate accelerates gradually（b）

2.6 开封河速率与温度的关系

河冰开河机制复杂，受热力、动力及河流形态（如多支流河流）等诸多因素影响，因此仅分析封河速率与温度的相关性。本文使用布尔津和哈巴河两个气象观测站的2001—2017年逐日平均气温，分析封河期平均温度与封河速率的相关性，结果显示，在封河期，布尔津站和哈巴河站观测到的平均温度随时间变化的斜率分别为0.155、0.124，均为正值具有上升的趋势［图9（a）、9（b）］，这与高沈瞳等［4］得到的额尔齐斯河流域年均温及夏、秋、冬温将继续保持升温的趋势的结果保持一致。图9（c）、9（d）显示，布尔津站封河期平均温度与封河速率的相关系数R=-0.179；哈巴河站封河期平均温度与封河速率的相关系数R=-0.299，封河期平均温度均与封河速率呈负相关，因此，随着封河期平均温度升高，封河速率变慢，本文认为封河期平均温度是封河速率变化的影响因子之一。

图9 两个气象观测站的封河期平均温度变化情况及封河期平均温度与封河速率的相关关系Fig.9 average temperature changes of two meteorological observation stations Burqin station（a）and Habahe station（b）during the river freeze-up period，and the correlation between the average temperature during the river freeze-up and the river freeze rate measured by Burqin station（c）and Habahe station（d）

3 结论与讨论

本文基于MOD09GQ数据，使用HANTS算法及TIMESAT软件，研究额尔齐斯河布尔津至斋桑泊河段的河冰物候，得到以下结论：

（1）通过HANTS算法和TIMESAT软件的处理，有效去除了云噪声对MODIS影像的影响，反演得到河冰物候信息。对比分析2015年2月2日（即最大雪/冰覆盖日）、3月15日、4月8日及4月12的反射率可知，随着河冰开河，河冰反射率明显降低，并且与南湾水文站观测数据的对比结果表明，MODIS数据所得与南湾水文站所得冰期持续时间的相关系数R=0.306呈正相关，本文研究结果与实测数据具有一致性，因此认为MOD09GQ的第二波段数据能有效监测河冰冰情变化。

（2）在额尔齐斯河布尔津至斋桑泊河段，开始封河时间和结束开河时间均有提前的趋势，最早开始封河时间为2016年9月29日，最晚开始封河时间为2013年11月13日，最早开始封河比最晚开始封河提前45天；最早结束开河时间为2013年3月28日，最晚结束开河时间为2002年4月28日，最早结束开河比最晚结束开河时间提前31天。封河与开河时间均提前，封河比开河提前的时间长，从而导致冰期持续时间有延长趋势。

（3）该河段封河速率随时间变化的斜率为-0.007，开河速率随时间变化的斜率为0.021，呈现出封河速率逐渐下降、开河速率逐渐加快的趋势。封河速率与布尔津站封河期平均温度的相关系数R=-0.179，与哈巴河站封河期平均温度的相关系数R=-0.299，封河速率与封河期平均温度均呈负相关，因此，封河期平均温度升高是引起封河期速率下降的一个因素。若未来该河段通航，则在封河期河流冻结速度变慢，可通航时间将延后，在开河期河冰加速消融，可通航时间将提前，每年总通航期延长，可带来更大经济效益。

本研究所得结果虽与实测数据所得结果相符，但尚有诸多不足之处，MODIS影像时间分辨率高，可有效监测河冰随时间的变化情况，但其空间分辨率较低，使用MODIS数据研究存在河宽不足250 m的河流，不能监测其具体结冰与破冰过程，并且可能存在混合像元的情况，为消除混合像元的影响，

使用整个河段的平均反射率进行研究，仅能描述所研究河段河冰整体随时间变化的趋势。此外，2015年4月16日高分1号影像缺失，并且直到4月24仅有一幅研究区上游影像，预览该影像显示研究区河段无河冰，未能观测到4月16日是否研究区全河段开河。额尔齐斯河属于跨国河流，可获取水文资料有限，仅能通过水文站监测所得冰期持续时间对本文结果进行验证。下一步可使用微波遥感数据监测河冰形成及破碎过程，结合使用无人机实际观察河流封河及开河过程，了解河冰具体变化过程，并验证MODIS数据所得结果。