GF-1遥感影像结合等高线消除云层干扰的连续水体重建法
2021-09-13徐国鑫刘林鑫1范亚洲1张企诺周佳奇1
张 珂，吴 南，徐国鑫，刘林鑫1，，范亚洲1，，张企诺，周佳奇1，，刘 挺
(1.河海大学水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；2.河海大学水文水资源学院，江苏 南京 210098； 3.长江保护与绿色发展研究院，江苏 南京 210098； 4.中国气象局-河海大学水文气象研究联合实验室，江苏 南京 210098； 5.陕西省引汉济渭工程建设有限公司，陕西 西安 710010)
随着空间技术的发展，通过遥感影像提取水体在洪涝灾害、水资源监测以及水环境保护等领域得到广泛应用[1-3]，水体提取已经成为遥感技术研究领域的重要分支之一[4]。目前常见的水体提取方法有单波段阈值法、谱间关系法、水体指数法、决策树法、密度分割法以及面向对象法等[5-9]，其中水体指数法是一种方便高效的水体识别法，应用也比较成熟[10]。由于卫星的分辨率不同，水体提取结果也有所差异，邓劲松等[11]基于SPOT高分辨率影像提出一种波段间的提取方法，有效地将居民地和水体分开；周成虎等[12]针对AVHRR影像提出了自动提取识别水体的描述模型，并在太湖、淮河等多个地区得到应用；王培培[13]利用ETM遥感影像数据，通过光谱特征分析和归一化植被指数识别水体，实现了水体的自动提取。以上研究所使用的数据源均来自国外遥感卫星，利用我国高分卫星进行水体提取的研究和算法还不多见。
高分一号卫星(GF-1)于2013年4月26日发射成功，是我国第一颗考核寿命大于5 a的卫星，卫星除装载有两台2 m分辨率全色和8 m分辨率多光谱相机外，还搭载了4台16 m分辨率多光谱相机。6台相机的配合使得我国首次实现了单颗卫星高分辨率与大幅宽技术的有效结合[14]，大大提高了我国高分辨率影像数据的自给率[15]，将为我国的科技发展带来不可估量的影响[16]。
受天气影响，遥感影像在地表水体提取过程中易受到云层的影响。遥感影像上受云层影响的区域覆盖了大面积的地表信息，使得图像信息不完善，给影像数据的处理和应用造成了极大的障碍[17-18]。如何消除遥感影像上云层的影响，至今仍是研究的热门话题。目前已有不少国内外相关专家和学者对此展开研究，研究结果表明：同态滤波法适用于大面积有薄云存在的遥感影像[19-21]，对厚云区域基本不产生效果；多源数据融合法使用不同传感器的遥感影像，将无云部分替代其他影像的云干扰部分[22]；时间序列傅立叶分析法简化了压缩时间序列遥感数据，对遥感影像中的云进行检测、去除和替代处理[23]。除此之外，将遥感影像与地形信息有效结合的研究也有很多[24-25]，王鑫蕊等[26]利用湖泊边界位于同一高度的特点，提出了一种利用湖泊的已知部分边界结合DEM数据，还原完整边界的自动方法，有效克服了云层的干扰。
屯溪流域为中小河流重点研发项目的典型示范流域，东方红水库是屯溪流域最大的一座中型水库，其库容、面积变化影响着该流域的应急防控技术。由于遥感数据周期间隔较长，无法实时获取水库的面积信息，且研究区内影像常年或多或少地受到云层的影响，降低了遥感数据的实用性。为了利用我国现有的光学遥感影像提高遥感数据的质量，本文以东方红水库为研究对象，提出一种GF-1遥感影像结合DEM消除云层干扰的连续水体重建法，对于连续水体局部受云层影响的影像，结合GF-1遥感影像和高精度等高线数据，通过地理配准还原遥感影像水体面积，从而实现云层干扰下水体的自动提取。
1 研究区概况与数据来源
1.1 研究区概况
东方红水库位于新安江水系二级支流东亭河(虞山溪)上(图1)，是以防洪、灌溉为主，结合发电、旅游供水等综合利用的中型水利枢纽工程。水库集水区域处在117°E～119°E、29°N～31°N之间，集水面积为60.0 km2，主河道长11.8 km。东方红水库总库容0.244 5亿m3，设计灌溉面积0.3万hm2，电站总装机容量400 kW，多年平均发电量131万kW·h。工程规模为安徽省重点中型水库，工程等别为Ⅲ等，拦河坝、溢洪道、引水涵洞，以及由导流隧洞改建的放水底孔等均为3级建筑物。
图1 研究区位置和地形Fig.1 Location and topography of the study area
1.2 数据来源
遥感数据来自中国资源卫星应用中心(http：//www.cresda.com/CN/index.shtml)的高分一号宽幅传感器(GF1_WFV)数据，影像空间分辨率为16 m。通过实地查勘的方式，收集了东方红水库高精度等高线数据及2015年水库详细运行资料。2015年为丰水年，东方红水库集水区域平均降水量为2 411.6 mm。通过2020年2月陆面过程资料分布式数据存档中心LP DAAC(Land Processes Distributed Active Archive Center)公布的1 rad/s分辨率NASADEM数据集产品(https：//lpdaac.usgs.gov/news/release-nasadem-data-products/)提取研究流域包含的水体区域。为体现水库在不同时期的蓄水面积信息，选取2015年非雨期(4、10月)、雨期(5、8、9月)共7景卫星影像。
2 研究方法
2.1 数据处理
2.1.1 遥感影像数据处理
选择中国资源卫星应用中心公布的1A级影像数据产品。为降低遥感影像在成像过程中大气传输、传感器本身、地球曲率等多种因素的影响，需要对原始遥感影像进行预处理。结合IDL(interactive data language)与ENVI(the environment for visualizing images)对数据进行批量预处理，预处理的方式主要包括：辐射定标、大气校正[27]、正射校正等。预处理结束后，将研究区域对应的矢量文件与遥感影像配准，完成研究区域的裁剪。2015年7景水体局部受云层影响的遥感影像数据见表1。
表1 高分一号卫星遥感影像数据集
2.1.2 高精度等高线数据处理
利用线性均分方式将东方红水库对应的高精度等高线数据进行插值，并选取多幅高保证率历史遥感影像，通过最大类间方差膨胀迭代法提取水体，从而得到水库不同水位值对应的水体边缘轮廓线，以此为参考，利用GIS软件对插值得到的高精度等高线数据进行人工解译，精准修正等高线数据，并确保高程数据完全覆盖水体范围。
2.2 水体提取与云层消除
2.2.1 最大类间方差膨胀迭代法提取水体
归一化差分水体指数(normalized difference water index，NDWI)法[9]利用遥感影像的绿波与近红外波两个特定波段进行归一化差值处理，从而突显影像中的水体。将经过预处理的遥感影像利用NDWI法初步提取水体。
最大类间方差法由日本学者Ostu等[28]提出，适用于双峰情况的阈值分割，可以实现自动求取阈值的目的。该方法根据图像的灰度值，将图像分为背景和目标两部分。分割阈值从最小灰度值到最大灰度值遍历每个值，不同的阈值分割使得背景和目标之间的方差不同，类间方差越大说明被阈值分割的两部分差别越大。若出现目标与背景的错分、漏分等现象，均会导致两部分差别变小。因此，利用类间方差最大时对应的阈值进行分割，被错分的概率应该是最小的。
利用NDWI法提取未受云层影响的局部水体，将初始阈值T0设为-0.02，若归一化差分水体指数INDWI
|Sk-Sk-1|<0.01 (k=1，2,3,…)
(1)
S0=Sk
(2)
式中：Sk——第k次迭代得到的水体面积；Sk-1——第k-1次迭代得到的水体面积；S0——最大类间方差膨胀迭代法提取到的最优水体面积。
2.2.2 基于高精度等高线数据的云层影响消除法
通过最大类间方差膨胀迭代得到东方红水库的局部水体范围，与处理后的高精度等高线数据进行地理坐标配准，可得到每个水体栅格对应的水位高程值。某一时刻连续水体水位应为唯一值，而由于云层的干扰以及误差影响导致遥感影像提取的边界水位并不唯一，故对云层干扰下由遥感数据得到的水体不同水位进行数学统计计算。
定义受云层影响遥感影像得到的局部水体为W类水体，局部水体边界网格以内的内部水体网格为W1，水库原始外轮廓边缘线对应的网格为W2，受云层影响而产生的边界线对应的水体网格为W3，其他非水体网格为W0。
若W2类水体网格对应的不同水位值共n个，对此进行均值计算：
(3)
SM⊕SV=SV
(4)
式中：SM——受云层影响遥感影像得到的局部水体面积；SV——还原水体水位值所包围的区域面积。两类水体经地理配准后可实现云层干扰下的水体还原提取。
3 结果与分析
在IDL + ENVI编程环境中实现遥感影像数据的预处理。将经过预处理的图像在MATLAB中进行最优阈值的确定，再通过GIS软件确定最优阈值对应的水体范围，与水体范围对应的高精度等高线数据进行地理配准，从而实现云层干扰下水体的还原处理。选取GF-1_2015-09-08景遥感影像(云量为11%)进行有云水体提取，各步骤处理结果如图2所示。
图2 水体GF-1遥感影像云消除过程Fig.2 Cloud removal process in water body extraction for the GF-1 satellite imageries
图2(a)为东方红水库集水区域2015-09-08GF-1_2015-09-08景遥感影像经过辐射定标、大气校正、正射校正、配准及研究区裁剪等预处理后的图像。由图2(a)可以看出，红色区域为水库水体区域，部分水体受到云层影响。图2(b)为云层干扰下通过最大类间方差膨胀迭代提取的东方红水库水体面积图，由于受云层影响较为严重，水库外轮廓线提取并不完整，且水库内部也缺失部分水体。图2(c)为东方红水库云干扰水体提取图与其高精度等高线数据的地理配准图，图中水体与等高线数据产生轻微偏移，水位值并不唯一。图2(d)为通过地理配准还原云层干扰下的东方红水库水体面积图，较为清晰地还原了水库的基本轮廓。
根据东方红水库289～297 m水位区间对应的水位-面积曲线(图3)以及观测资料记录的每日水库实测水位值(表2)，可以推导出对应时间水库的水体面积，以此作为SV的参考值。
图3 东方红水库水位-面积关系曲线Fig.3 Water level-area relationship curve of Dongfanghong Reservoir
从表2可以看出，东方红水库GF-1_2015-09-08景遥感影像实测水位值h为291.35 m，根据水位-面积关系曲线得到的参考面积S为1.995 6 km2。然而由于云层的干扰，根据最大类间方差膨胀迭代得到的水体面积为1.060 7 km2，面积差值为0.934 9 km2，差距较大。按照本文提出的通过地理配准还原云层干扰水库水体面积，可得到水库水位还原值为291. 95 m，还原后库区水体面积为2.042 km2，面积差值为0.046 4 km2，可以看出，该方法简单有效地还原了水库水体面积。
表2 消除云层干扰后的库区水体遥感提取面积
另选取东方红水库2015年6景影像并使用上述方法进行云消除处理。图4为最大类间方差膨胀迭代提取的水体范围与其等高线数据的地理配准图，蓝色区域为S0。从图4可以看出云层干扰对东方红水库水体提取范围影响较大，严重时水库外轮廓基本无法确定。除此之外还可以看出，库区水体外轮廓线对应的水位值并不唯一，因此需要通过数学统计计算方式确定需要还原的水位值。
图4 地理配准后的卫星影像Fig.4 Satellite imageries after geographic referencing
图5为经地理配准后还原的水体，与等高线数据相结合经均值计算后，水库对应的水位值可唯一确定，最终结果表明该方法简单有效地消除了云层的干扰，较准确地还原了东方红水库在非雨期(4、10月)、雨期(5、8、9月)云层影响下的水体。从图5(a)(b)(f)可以看出，由于水库在4、10月来水量较少，随着下游水量的补给，库区水体范围逐渐减小，提取的水体范围较小；图5(c)(d)(e)表明，5、8、9月降水量较多，水库来水量大于出水量，蓄水量逐渐增多，故水体面积增大，水位升高，库区之间连接的细小河道蓄水位也逐渐升高。
图5 水体还原结果Fig.5 Results of extracted water bodies
从表2可以看出，S0与S相比差距较大，其中GF-1_20150426景遥感影像水体面积相差最大，达1.368 4 km2，难以反映真实水体。选取的7景影像经地理配准后得到的水体还原值与实测水位参考值接近，还原后的水体水位平均差值为0.438 6 m，有效水体面积平均差值为0.080 9 km2，平均相对误差在5%以内，较为精确地还原了库区的水体。其中GF-1_20150512景遥感影像水体库区面积差值仅为0.008 3 km2，较好地消除了云层干扰，还原了水体水位值及对应的面积。
从表2还可以看出，GF-1_20150401、GF-1_20150426和GF-1_20151003景遥感影像还原的水位值与实测水位值差距略微偏高，这可能是由于所用影像分辨率较低，为16 m，在预处理过程中水体范围整体产生轻微偏移，使得提取的水体存在部分误差，且该方法对地形数据精度要求较高，本文将1 m等高线经线性插值、人工修订后得到的等高线与提取的遥感信息进行地理配准，不同数据源之间易出现误差叠加现象。
4 结 语
高分一号卫星遥感影像受云层影响限制了数据的实用性，且在云层干扰下通过最大类间方差膨胀迭代法得到的水体并不完整，面积平均差值为1.009 4 km2，无法准确反应水体的真实信息，难以直接用于陆地表面水体的提取。通过本文提出的地理配准还原云层干扰下水体的方法，还原后的水体水位平均差值为0.438 6m，有效水体面积平均差值为0.080 9 km2，最小仅为0.008 3 km2，平均相对误差在5%以内，有效地克服了云层干扰，较为精确地还原了库区水体。该方法简单有效，与高精度等高线数据相结合，通过简单的图像计算便可克服云层的干扰，提高了卫星遥感影像的实用性，可广泛用于地表水资源的动态提取，为水资源调查及地表水体动态监测等领域提供高质量的遥感影像数据。
通过本文方法得到的水体面积虽然有效克服了云层的干扰，但该方法对地形数据精度要求较高，对于缺少水下高程数据的水体来说，水位动态变化值常在某范围内上下变动，故可根据多景无云历史遥感影像(包含丰水期、枯水期较多高程值变化)得到的水位值对水下空缺的高程值进行填补，再进一步进行人工修正，但精度方面问题有待提高。除此之外，对于一些复杂的水体情况，如湖心岛等，该方法需要结合人工解译才能进行水体提取，对此仍需要进一步的深入研究。
