吴 杰,董小涛,张 珂,3,4,李 曦,吴 南

(1.河海大学水文水资源学院,江苏 南京 210098; 2.水利部综合事业局,北京 100032;3.长江保护与绿色发展研究院,江苏 南京 210098; 4.中国气象局水文气象重点开放实验室,江苏 南京 210098)

消落带作为一种特殊的湿地生态系统,指由于季节或人为控制因素导致水位消涨而在水库周边形成的陆地和水生生态系统功能相互作用的巨大隔离带,是库区生态保护的重要对象[1-3]。水位波动是对湖泊生态系统影响最为重要的因素[4],也是形成消落带的主要原因。与自然形成的河岸带不同,水库消落带的水位变化主要受水库实际调度影响,水库调蓄造成的库区水位周期性波动变化会对消落带的生态结构造成重大影响。目前对消落区的研究主要集中在土地利用变化[5-6]、消落带分区等方面,对消落带时空变化特征分析的研究较少。

淹没频率作为一项重要的水文特征,反映了空间位置上水体淹没时间的长短及次数,是影响湿地生态系统的重要因素[7],同时会影响消落带的生态结构和生态类型划分。张丽丽等[8]通过计算鄱阳湖自然保护区湿地的水深和淹水频率,分析了植被群落对这2种水文特征的耐受性和敏感性。许超等[9]基于MODIS遥感数据计算淹没频率,分析其变化;谷娟等[10]基于混合像元分解模型分析了鄱阳湖的淹没频率变化,探究时间序列上湿地植被特征的空间响应。水位的季节性波动变化导致水库周边形成不同的淹没频率,湿地中的植被也受到不同程度的水淹影响,对水位变化不适应的植被大量消亡,造成该湿地区域植被群落减少,环境的稳定性受到破坏,水土流失、水质恶化等损害生态系统的问题层出不穷。因此,综合分析消落带淹水范围及淹水频率的时空变化特征,加强对消落带水体分布的高精度、高频次监测具有重要意义。

遥感信息技术作为一种可以高效获取、分析及处理空间地理信息的技术手段,已广泛应用于水文领域的趋势变化分析及宏观监测[11-15]等领域。高分一号卫星由我国自主研发并于2013年4月发射成功,有效载荷为2台2m分辨率全色、1台8m分辨率多光谱相机(PMS)和4台16m分辨率多光谱宽幅相机(WFV),具有高空间分辨率、多相机等特点,满足了我国对高时空分辨率遥感数据的应用需求[16]。目前常见的遥感影像水体信息提取方法包括单波段阈值法[17]、多波段谱间关系法[18]、水体指数法及分类法[19- 20]等。

为研究消落带生态类型划分,本文提出一种线性插值模拟计算水体淹没频率的方法,并结合遥感影像水体信息提取[21-23]和几何校正[24-25]技术,分析了三河口水库蓄水阶段的水体时空变化特征,并根据模拟水体的精度反映淹没频率结果的可靠性,以期提高淹没频率计算的可靠性。

1 研究方法

1.1 遥感数据处理及精细配准

本文采用的遥感数据为中国资源卫星网站发布的高分一号卫星WFV16 m分辨率遥感影像。在进行后续数据分析及应用前,利用ENVI二次开发工具IDL批量预处理不同日期遥感影像,步骤包括辐射定标、大气校正[26]、正射校正等。为了消除因传感器自身因素、大气条件、太阳位置和角度等引起的遥感影像空间位置误差,选取固定的地面建筑物作为参照物进行精细配准,使配准后的影像空间位置差异在误差允许范围内[24]。根据研究区矢量文件对预处理后的遥感影像进行裁剪得到用于水体提取的源数据。

1.2 最大类间方差迭代水体提取法

最大类间方差法又称为大津法(Otsu),由学者Otsu[27]于1979年提出,是一种简单可行、应用广泛的图像分割阈值算法。该算法可以找到某一灰度值,即阈值,将图像分为前景和背景。本文根据归一化差分水体指数(normalized difference water index,NDWI)将图像分为水体和非水体。阈值为使前景和背景的类间方差最大的灰度值,类间方差越大说明前景与背景灰度值的离散程度越大,同时也说明类内方差越小,即前景和背景内部灰度值离散程度越小,故此时的阈值选取也更加合理。

基于改进的最大类间方差迭代水体提取法[11],通过迭代不同膨胀算子确定最优提取结果。首先通过经验阈值和边缘检测法初步筛选水体信息并利用连通域标记法去除其他无关水体,然后迭代多个膨胀算子并基于形态学膨胀法建立库区水体的缓冲区,再对缓冲区内区域应用最大类间方差法,根据判断准则自适应地确定水体与非水体的分割阈值。自适应最优阈值选择公式为

Smin=min{|Si-Si-1|} (i=1,2,3)

(1)

式中:Smin为最小的水体提取面积差值;Si、Si-1分别为第i次、第i-1次迭代得到的水体面积。

1.3 连通域标记法

连通域指的是图像中由相同并且相邻的像素值组成的图像区域。不同的相邻像素点形成不同的区域,这些不同区域代表了不同的连通区域,像素的邻接关系决定了像素与周边的连通关系,通常分为四邻接和八邻接关系[28]。本文图像为二值图像,最小单位为像素,采用四邻接关系,即考虑像素的上、下、左、右等4个方向的连接状态,通过遍历图像,记下每一行或列中团和标记的等价对,并通过等价对对原来的图像进行重新标记,步骤分为横向合并、纵向合并和逆向合并。

通过连通域区域标记算法对各连通水体进行标记,对研究区内不同水体进行区分,选择最大连通域作为库区水体范围,可以有效去除图像中存在的其他水体和噪音。

1.4 水库淹没频率计算方法

水库淹没频率是指某个位置上被水淹的时间与总研究时间的比值。受遥感影像时间分辨率、云体干扰、数据质量等因素干扰,无法获得日尺度的水体时空变化特征。本文选择影像在时间尺度上保证相邻两次日期高分一号卫星数据的水库水体时空变化呈单调变化趋势,对缺少数据的日期,基于线性插值模拟结果构建淹没频率计算公式:

(2)

式中:P为对应像元在研究时间内的淹没频率;Nj为第j+1期影像该位置的水体划分,水体与非水体分别对应为1、0;ΔT为两期影像间隔时间;T为总研究时间。

2 实例验证

2.1 工程概况

引汉济渭工程为陕西省内重大跨流域调水工程,三河口水利枢纽是引汉济渭工程的调蓄中枢,目前已完成水库初期蓄水阶段,消落带正在逐渐形成。三河口水利枢纽位于佛坪县与宁陕县境交界、汉江一级支流子午河中游峡谷段,其坝址位于大河坝乡三河口村下游约2km处。三河口水利枢纽主要由拦河大坝、泄洪放空系统、供水系统和连接洞等水工建筑物组成。

三河口水库回水淹没范围涉及陕西省宁陕县,库区上游左岸有椒溪河、蒲河和汶水河等支流汇入。该水库库容较小,是整个引汉济渭工程中具有较大水量调节能力的核心项目,承担供水、调蓄、发电以及满足下游综合用水要求等多项任务。

根据三河口水库的蓄泄运行规律,水库水位在正常蓄水位643m至库正常运行死水位558.0m之间变动。图1展示了三河口水库最大淹没范围1000m缓冲区。

图1 研究区示意图Fig.1 Schematic diagram of the study area

2.2 结果与分析

2.2.1 库区水位变化趋势

选取2020年1月22日到2021年11月8日共25幅预处理后的高分一号卫星遥感影像,应用Otsu法选取阈值得到水体提取结果。自2020年1月22日至2020年2月29日水库始终处于最低水位,2021年9月8日后达到最大水位,涨水阶段共历时约18月,预计消落带总面积为13.95km2。对遥感影像进行目视解译,得到三河口水库水体面积的变化阶段:2020年2月以前,水库维持在死水位不变;2020年2月至2021年2月,水库水位第一次上涨后维持固定水位,涨落周期为11月;2021年2—6月,水库水位在大幅上涨后维持固定水位,涨落周期为4月;2021年6—10月水库水位继续上涨至正常蓄水位,之后水位维持在正常蓄水位。

图2为时间序列上Otsu法提取水体的变化趋势,线性拟合R2达到0.58468,总体呈现稳步上升趋势。20201129景影像水体提取结果较相邻日期波动较大,结合目视解译遥感影像发现将大坝周围建筑及山体阴影误分为水体导致水体提取结果较实际结果偏大;20210117、20210605景影像水体提取结果较周围日期向下波动则是由于受云层影响水体上游末端细小水体未充分提取,导致水体提取结果较实际结果偏小。三河口水库周边存在部分建筑和山体阴影,而NDWI指数对山体阴影、建筑物和水体的区分能力较弱,存在将其误分为水体导致部分日期水体面积偏大的情况,呈现为图2中蓄水过程线的波动。

图2 三河口库区水体面积时间序列Fig.2 Time series of water body area in Sanhekou Reservoir

由图3可知,水体边界整体沿水系扩展,水体面积增加,水位由558m上涨至643m,落差达85m,符合水库实际调度情况。其中,20210219景影像之前库区水体面积变化较小,存在将坝址、山体阴影错分的现象;20210219景影像后随着水库大量蓄水,水位大幅上升,水体面积增加,水体与周边噪音的区分程度较好,对水体的误提取情况明显减少。

图3 2020—2022年三河口库区水体提取结果Fig.3 Results of extracted water body area of Sanhekou Reservoir from 2020 to 2022

2.2.2 水库空间淹没频率特征

对提取的2020—2022年三河口水库水体边界,在进行大气校正等预处理手段后,通过精细几何校正使影像之间地物的误差控制在误差允许范围内。将处理后的影像进行叠置分析并按时间序列进行插值处理,计算年淹没频率。

图4(a)为得到的空间上16m分辨率每个像元的淹没频率,反映了库区空间淹没状态,淹没程度沿三河口水库三条河流上游向库区逐渐增加,较真实地反映了库区淹没状态。但永久淹没区较死水位对应的水面面积偏小,分析误差原因主要为不同遥感影像间空间位置差异。通过人工修正按淹没频率(0, 33.3%]、(33.3%, 66.6%]、(66.6%, 100%)、100%将三河口水库水位变动区域划分为4个淹没等级(图4(b)):经常出露、半淹半露、经常淹没和永久淹没。可以看出坝址周围建筑及山体阴影的误分影响了淹没频率的计算。

图4 三河口水库1000m缓冲区内的淹没频率特征及淹没频率分区Fig.4 Characteristics and zoning of inundation frequency with 1000m buffer zone in Sanhekou Reservoir

从表1可以看出淹没频率面积占比随淹没程度的增加逐渐递减,空间上表现为由内向外递减的趋势。由于研究期限内水库正处于初期蓄水阶段,淹没频率区间为(0, 33.3%]的区域占比最大,面积达到11.72km2,占淹没区面积的81.84%,是淹没区的主体部分;(33.3%, 66.6%]、(66.6%, 100%)、100%的区域面积分别为1.55、0.68、0.37km2,占淹没区面积的10.82%、4.75%、2.58%。

表1 缓冲区淹没频率统计

2.2.3 水库淹没频率可靠性分析

选择存在遥感影像的日期进行验证,基于淹没频率空间特征结果,模拟得到水体分布结果。淹没频率结果的可靠性通过模拟得到的水体分布结果的精度来反映。用于验证的遥感影像需要满足未参与频率计算且处于水库调蓄阶段内。

由图5可知,20210312与20210323、20210708、20210717水体提取结果基本一致,说明它们处于水库上涨的同一阶段。结合人工目视解译实际遥感影像,分别选取100个水体和非水体像元构建混淆矩阵计算总体精度和Kappa系数,结果见表2。

表2 水体模拟结果精度

图5 不同日期模拟水体分布Fig.5 Simulated water body distribution

由表2可知, 6次水体模拟结果的总体精度和Kappa系数均值分别达到86.5%和73%,结果较为满意,说明由本文方法得到的淹没频率特征符合该地区的水库调蓄情况,具有较强的可靠性。

3 结 语

本文采用线性插值法计算了不同时段的水淹状态,模拟了蓄水时段的水体变化,结合高分辨率遥感数据,得到库区淹没频率的时空分布特征。根据淹没程度不同分成经常出露、半淹半露、经常淹没、长久淹没4种淹没区,空间上呈现为由内向外递减,面积占比分别为81.84%、10.82%、4.75%、2.58%,在水库进入正常调蓄阶段后,为消落带形成后的生态类型划分提供了借鉴思路。本文选择未参与计算且存在遥感数据的日期对水体模拟结果进行精度验证,总体精度最大值为90%,最小值为82%,均值为86.5%,Kappa系数最大值为80%,最小值为64%,均值为73%,说明淹没频率时空分布特征可靠性较好。总体上该频率计算方法简单有效,在水库调蓄情况下较为适用,精细化了淹没频率的计算。但强降雨等自然因素的干扰会影响其适用性;水体提取结果的准确性也会影响频率计算的结果;在小型水库由于对几何校正的要求变高也会导致精度的降低,对此在未来的研究中仍需进一步深入探讨。