CN103363962B - 一种基于多光谱影像的湖泊水储量遥感估算方法 - Google Patents

Abstract

本发明提供一种基于多光谱影像的湖泊水储量遥感估算方法，该方法包括：根据湖泊水体在多光谱遥感影像不同波段上的光谱响应特性提取反映湖泊水面分布特征的水体指数；根据提取的水体指数获取湖泊水体边界，并将实测的湖泊水位信息赋给湖泊水体边界，根据不同时期的湖泊水体边界生成湖泊等水位线数据集；根据湖泊等水位线数据集模拟湖泊水下地形；根据模拟的湖泊水下地形和实测的水位数据计算不同时期的湖泊水储量。该方法利用多光谱卫星遥感数据对湖泊水下地形进行监测，由于使用的卫星遥感数据基本上可以免费获取，而且遥感数据空间覆盖范围广、更新速度快，因此，相比传统的野外实测方法，本发明实施例提出的方法具有监测成本低，更新方便及能够在大范围内进行推广应用的优点。

Description

一种基于多光谱影像的湖泊水储量遥感估算方法

技术领域

本发明涉及卫星遥感和对地观测领域，尤其涉及一种基于多光谱影像的湖泊水储量遥感估算方法。

背景技术

湖泊（包括天然湖泊和人工库塘）是地表水资源的最主要表现形式。除冰川和永久雪盖外，湖泊水是第二大地表水资源类型，它是影响全球海平面变化的主要因素之一。湖泊演化与环境变化密切相关：一方面，湖面波动引起的水量变化和沉积环境改变可敏感地指示气候变化；另一方面，湖泊水面积的增减会改变下垫面条件，从而对气候变化产生影响。然而长期以来，受调查和观测资料、模拟和分析方法的限制，人们无法精确获取大范围湖泊水资源储量并描述其动态变化特征，严重影响了人们对陆表水储量变化与全球或区域水循环加剧、全球海平面升降的关系等科学问题的认识。在水资源管理应用方面，这一现状也限制了人们对近年来全球范围内频发的旱涝灾害事件的应对时效。

因此，需要采用有效的方法对湖泊的水储量进行准确估算。现有的估算方法是：以实测水下地形数据为基础构建的水位-容积曲线法。该方法主要是根据实测水下地形数据和水位测量数据采用相应的数学模型进行体积的计算，从而获得水位-容积曲线方程，在获知湖泊实时水位后，就能根据水位-容积曲线得知该水位相对应的容积值。

上述方法的应用受限于库容计算数学模型、水位测量空间代表性和水下地形这三个方面资料的获取。

常用的库容计算方法有断面法、等高线容积法、方格网法和三角网格法，受湖泊形态、大小和复杂程度的影响，不合适的数学模型会带来较大的计算误差。这一问题可以利用当前的高性能计算资源，通 过增加模型微分量级加以解决。对于水位观测站测量结果空间代表性差的问题，则可以通过优化观测站点位置和增加观测站点数量来解决。

而对于水下地形资料来说，一般地区都有早期l:5000或1:10000地形图，可通过收集扫描矢量化，采集高程数据并生成数字地形模型（DTM）。此方法生成的DTM数据工期短，但这些基础资料的收集一般较为困难，获取的库容数据精度易受地形图本身精度的影响，且无法反映制图后期上游水沙及人类活动对水下地形改造的影响。解决这一问题的办法是综合利用3S技术，即GPS、GIS、RS，定期开展湖泊水下地形测量，如综合利用GPS和激光雷达数据的水下地形制图，基于GPS野外调查数据的水下地形模拟，以及直接利用激光雷达获取的数字高程模型（DEM）模拟湖泊水下地形等。这些方法虽然能部分解决水下地形数据获取及更新问题，但受测量周期和费用的限制，无法进行大范围推广应用。

发明内容

（一）要解决的技术问题

本发明的目的是提出一种湖泊水储量估算方法，使得该方法能够在大范围内进行推广应用。

（二）技术方案

为了解决上述技术问题，本发明提出了一种基于多光谱影像的湖泊水储量遥感估算方法，该方法包括：

S1、根据湖泊水体在多光谱遥感影像不同波段上的光谱响应特性提取反映湖泊水面分布特征的水体指数；

S2、根据提取的水体指数获取湖泊水体边界，并将实测的湖泊水位信息赋给湖泊水体边界，根据不同时期的湖泊水体边界生成湖泊等水位线数据集；

S3、根据湖泊等水位线数据集模拟湖泊水下地形；

S4、根据模拟的湖泊水下地形和实测的水位数据计算不同时期的 湖泊水储量。

其中，在步骤S1前，首先进行遥感数据的收集和遴选；并将收集和遴选的遥感影像进行辐射纠正和几何纠正；然后根据纠正后的遥感影像提取反映湖泊水面分布特征的水体指数。

其中，水体指数的提取可以采用两种途径实现：

根据湖泊水体在Landsat MSS或HJ-1A/B影像不同波段上的光谱响应特征提取归一化差异水体指数；

或者，

根据湖泊水体在Landsat TM/ETM+影像不同波段上的光谱响应特征提取增强型归一化差异水体指数。

其中，归一化差异水体指数NDWI采用下式计算：

$\mathrm{NDWI}=\frac{{\mathrm{\ρ}}_{\mathrm{Green}}-{\mathrm{\ρ}}_{\mathrm{NIR}}}{{\mathrm{\ρ}}_{\mathrm{Green}}+{\mathrm{\ρ}}_{\mathrm{NIR}}}$

其中，ρ_Green为Landsat MSS或HJ-1A/B影像上的绿光波段，ρ_NIR为Landsat MSS或HJ-1A/B影像上的近红外波段。

其中，增强型归一化差异水体指数MNDWI采用下式计算：

$\mathrm{MNDWI}=\frac{{\mathrm{\ρ}}_{\mathrm{Green}}-{\mathrm{\ρ}}_{\mathrm{SWIR}}}{{\mathrm{\ρ}}_{\mathrm{Green}}+{\mathrm{\ρ}}_{\mathrm{SWIR}}}$

其中，ρ_Green为Landsat TM/ETM+影像上的绿光波段，ρ_SWIR为Landsat TM/ETM+影像上的短波红外波段。

优选地，根据提取的水体指数采用改进的最大类间方差阈值法提取湖泊水面栅格数据；通过ArcMap9.3软件中栅格转矢量工具将提取的湖泊水面栅格数据转换成湖泊水体边界矢量数据。

其中，根据提取的水体指数C的灰度范围{0,1,…,d-1}，将C分割为C₁和C₂两类，分割的灰度阈值为t，即湖泊的水体边界阈值，t满足：

$t=\arg \underset{0\≤t\≤d-1}{\max }\frac{S^2}{S_1^2+S_2^2}$

其中，S²为C₁和C₂的类间方差，S₁ ²为C₁的类内方差，S₂ ²为C₂的类内方差；

其中，

S²=P₁(A₁-A)²+P₂(A₂-A)²

$S_1^2=\frac{{\mathrm{\Σ}}_{i=0}^t{(i-A_1)}^2{p}_i}{P_1}$

$S_2^2=\frac{{\mathrm{\Σ}}_{i=t+1}^{d-1}{(i-A_2)}^2{p}_i}{P_2}$

其中，P₁为C₁中像元数占水体指数C中总像元数的比例，P₂为C₂中像元数占水体指数C中总像元数的比例，i为像元灰度值，A为C的像元灰度平均值，A₁为C₁的像元灰度平均值，A₂为C₂的像元灰度平均值，p_i为灰度值为i的像元数占水体指数C中总像元数的比例。

其中，根据湖泊等水位线数据集采用不规则三角网模拟湖泊水下地形，即湖泊水下地形由一系列互不重叠的不规则三角形组成的网状结构来表示，每个三角形的每个节点包括独立的经纬度坐标、水深值和坡度角；

根据模拟的水下地形和实测的水位数据将湖泊水体简化成一系列三棱柱的集合，整个湖泊水体的体积则为各三棱柱体积的和，计算公式为：

$V=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{S}_i\frac{h_i+h_{i+1}+h_{i+2}}{3}$

其中，n为三棱柱的个数，S_i为第i个三棱柱水体上表面的面积，h_i、h_i+1和h_i+2为三棱柱三条侧棱所在位置的水深值。

优选地，根据湖泊等水位数据集采用ArcMap9.3软件的三维分析工具中的TIN创建模块模拟湖泊水下地形。

优选地，湖泊水储量利用ArcMap9.3软件三维分析工具中的面积-体积统计模块进行计算。

（三）有益效果

本发明的实施例中利用多光谱卫星遥感数据对湖泊水下地形进行监测，由于使用的卫星遥感数据基本上可以免费获取，而且遥感数据的空间覆盖范围广、更新速度快，因此，相比传统的野外实测方法，本发明实施例提出的方法具有监测成本低，更新方便及能够在大范围内进行推广应用的优点。

另外，在湖泊水下地形遥感监测结果为基础的湖泊水储量遥感估算方法，能有效的消除不同时期湖泊水下地形本身变化对估算结果的影响，其精度相较传统方法更为可靠。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的实施例中基于多光谱遥感影像的湖泊水储量估算方法概念图；

图2为本发明的实施例中基于多光谱遥感影像的湖泊水储量估算方法的流程图。

具体实施方式

下面结合附图和实施例对本发明的实施方式作进一步详细描述。以下实施例用于说明本发明，但不能用来限制本发明的范围。

为了能够更好地掌握大范围内湖泊水储量及其变化情况，本发明实施例提出了一种基于多光谱影像的湖泊水储量遥感估算方法，参见图1，该方法包括：

S1、根据湖泊水体在多光谱遥感影像不同波段上的光谱响应特性提取反映湖泊水面分布特征的水体指数；

S2、根据提取的水体指数获取湖泊水体边界，并将实测的湖泊水 位信息赋给湖泊水体边界，根据不同时期的湖泊水体边界生成湖泊等水位线数据集；

S3、根据湖泊等水位线数据集模拟湖泊水下地形；

S4、根据模拟的湖泊水下地形和实测的水位数据计算不同时期的湖泊水储量。

本发明的上述实施例中，利用多光谱卫星遥感数据对湖泊水下地形进行监测，由于使用的卫星遥感数据基本上可以免费获取，而且遥感数据的空间覆盖范围广、更新速度快，因此，相比传统的野外实测方法，本发明实施例提出的方法具有监测成本低，更新方便及能够在大范围内进行推广应用的优点。

本发明的上述实施例中，由于不同卫星可以提供不同时期、不同质量的多光谱遥感影像，因此，方法使用前需要对所需要的遥感影像进行收集和遴选；另外，为了消除云雾覆盖、卫星姿态变化和地形起伏等对卫星成像的影响，在使用遥感影像之前需要对遥感影像进行辐射纠正和几何纠正，然后根据纠正后的遥感影像提取可以反映湖泊水面分布特征的水体指数，参见图2。

本发明的上述实施例中使用的水体指数包括归一化差异水体指数NDWI和增强型归一化差异水体指数MNDWI，其中NDWI适用于Landsat MSS和HJ-1A/B影像，计算公式如下：

$\mathrm{NDWI}=\frac{{\mathrm{\ρ}}_{\mathrm{Green}}-{\mathrm{\ρ}}_{\mathrm{NIR}}}{{\mathrm{\ρ}}_{\mathrm{Green}}+{\mathrm{\ρ}}_{\mathrm{NIR}}}---\left(1\right)$

式中，ρ_Green为Landsat MSS或HJ-1A/B影像上的绿光波段，ρ_NIR为Landsat MSS或HJ-1A/B影像上的近红外波段。

而MNDWI适用于Landsat TM/ETM+影像，其计算公式如下：

$\mathrm{MNDWI}=\frac{{\mathrm{\ρ}}_{\mathrm{Green}}-{\mathrm{\ρ}}_{\mathrm{SWIR}}}{{\mathrm{\ρ}}_{\mathrm{Green}}+{\mathrm{\ρ}}_{\mathrm{SWIR}}}---\left(2\right)$

式中，ρ_Green为Landsat TM/ETM+影像上的绿光波段，ρ_SWIR为Landsat TM/ETM+影像上的短波红外波段。

本发明的上述实施例中，根据提取的水体指数采用改进的最大类间方差阈值法提取湖泊水面栅格数据；通过ArcMap9.3软件中栅格转矢量工具将提取的湖泊水面栅格数据转换成湖泊水体边界矢量数据。

本发明的上述实施例中，改进的最大类间方差阈值法的基本原理如下：

根据提取的水体指数C的灰度范围{0,1,…,d-1}，将C分割为C₁和C₂两类，分割的灰度阈值为t，即湖泊的水体边界阈值，t满足：

$t=\arg \underset{0\≤t\≤d-1}{\max }\frac{S^2}{S_1^2+S_2^2}---\left(3\right)$

其中，S²为C₁和C₂的类间方差，S₁ ²为C₁的类内方差，S₂ ²为C₂的类内方差；

其中，

S²=P₁(A₁-A)²+P₂(A₂-A)²           (4) 

$S_1^2=\frac{{\mathrm{\Σ}}_{i=0}^t{(i-A_1)}^2{p}_i}{P_1}---\left(5\right)$

$S_2^2=\frac{{\mathrm{\Σ}}_{i=t+1}^{d-1}{(i-A_2)}^2{p}_i}{P_2}---\left(6\right)$

其中，P₁为C₁中像元数占水体指数C中总像元数的比例，P₂为C₂中像元数占水体指数C中总像元数的比例，i为像元灰度值，A为C的像元灰度平均值，A₁为C₁的像元灰度平均值，A₂为C₂的像元灰度平均值，p_i为灰度值为i的像元数占水体指数C中总像元数的比例。

本发明的上述实施例中，根据遥感影像获取湖泊等水位数据集后，当湖泊水下地形比较复杂时，可采用不规则三角网（TIN）的方法模拟湖泊水下地形，即湖泊水下地形由一系列互不重叠的不规则三角形组成的网状结构来表示，每个三角形的每个节点包括独立的经纬度坐标、水深值和坡度角；水下地形的模拟过程可以采用ArcMap9.3软件的三维分析工具中的TIN创建模块进行模拟，将得到的湖泊等水位数据集作为软件的输入数据，运行软件后，得到使用不规则三角网模 型模拟的湖泊水下地形。

以得到的湖泊水下地形数据为基础，结合实测的水位数据可将湖泊水体简化成一系列三棱柱的集合，整个湖泊水体的体积则为各三棱柱体积的和，计算公式为：

$V=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{S}_i\frac{h_i+h_{i+1}+h_{i+2}}{3}---\left(7\right)$

其中，n为三棱柱的个数，S_i为第i个三棱柱水体上表面的面积，h_i、h_i+1和h_i+2为三棱柱三条侧棱所在位置的水深值。

上述对湖泊水储量的估算可采用ArcMap9.3软件三维分析工具中的面积-体积统计模块进行计算，通过输入获得的湖泊水下地形和湖泊水位数据，运行软件后，即可得到湖泊水储量值。

当湖泊水下地形相对简单时，可采用方格网模拟湖泊水下地形。

因此，当需要对某个湖泊的水储量进行估算时，可按照下述方法进行：

1、获取覆盖湖泊区域的不同时间段的多光谱卫星遥感影像，如Landsat MSS/TM/ETM+和HJ-1A/B；

2、对获取的遥感影像进行辐射校正和几何校正；

3、利用校正后的影像数据计算湖泊区不同时期的水体指数，包括NDWI和MNDWI，采用（1）和（2）式进行计算；

4、采用改进的最大类间方差阈值法从水体指数数据集中提取湖泊水面栅格数据集，再利用ArcMap9.3软件中的栅格转矢量工具将湖泊水面栅格数据集转换成湖泊水体边界矢量数据集，再将收集的不同时期的湖泊水位数据赋给对应时间段的水体边界矢量，即可获得不同时期的湖泊区等水位线数据集；

5、根据获得的等水位线数据集，采用不规则三角网模拟方法，利用ArcMap9.3软件三维分析工具中的TIN创建模块模拟水下地形；

6、根据模拟的水下地形和实测的水位数据，利用ArcMap9.3软件三维分析工具中的面积-体积统计模块计算湖泊的水储量。

本发明的上述实施例中，对于不同的区域，采用的卫星遥感影像数据可能不同，如CBERS、ASTER、MODIS、SPOT等遥感影像；而实测湖泊水位数据的使用方法也不一样，本发明的实施例中使用的是湖泊区实测水位的均值，而对于水位观测数据丰富的区域，同一水体边界上可以采用不同水位值。

因此，本发明的有益效果是：

本发明的实施例中利用多光谱卫星遥感数据对湖泊水下地形进行监测，由于使用的卫星遥感数据基本上可以免费获取，而且遥感数据空间覆盖范围广、更新速度快，因此，相比传统的野外实测方法，本发明实施例提出的方法具有监测成本低，更新方便及能够在大范围内进行推广应用的优点。

另外，以湖泊水下地形遥感监测结果为基础的湖泊水储量遥感估算方法，能有效的消除不同时期湖泊水下地形本身变化对估算结果的影响，其精度相较传统方法更为可靠。

本发明的实施例是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显而易见的。选择和描述实施例是为了更好说明本发明的原理和实际应用，并且本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。

Claims (8)

1.一种基于多光谱影像的湖泊水储量遥感估算方法，其特征在于，包括：

S1、根据湖泊水体在多光谱遥感影像不同波段上的光谱响应特性提取反映湖泊水面分布特征的水体指数；

S2、根据提取的水体指数获取湖泊水体边界，并将实测的湖泊水位信息赋给湖泊水体边界，根据不同时期的湖泊水体边界生成湖泊等水位线数据集；

S3、根据湖泊等水位线数据集模拟湖泊水下地形；

S4、根据模拟的湖泊水下地形和实测的水位数据计算不同时期的湖泊水储量；

步骤S2中根据提取的水体指数获取湖泊水体边界包括：

根据提取的水体指数采用改进的最大类间方差阈值法提取湖泊水面栅格数据；

通过ArcMap 9.3软件中栅格转矢量工具将提取的湖泊水面栅格数据转换成湖泊水体边界矢量数据；

其中，所述根据提取的水体指数采用改进的最大类间方差阈值法提取湖泊水面栅格数据包括：

根据提取的水体指数C的灰度范围{0,1,…,d-1}，将C分割为C₁和C₂两类，分割的灰度阈值为t，即湖泊的水体边界阈值，t满足：

$t=\arg \underset{0\≤t\≤d-1}{\max }\frac{S^2}{S_1^2+S_2^2}$

S²为C₁和C₂的类间方差，S₁ ²为C₁的类内方差，S₂ ²为C₂的类内方差；

S²＝P₁(A₁-A)²+P₂(A₂-A)²

$S_1^2=\frac{{\mathrm{\Σ}}_{i=0}^t{(i-A_1)}^2{p}_i}{P_1}$

$S_2^2=\frac{{\mathrm{\Σ}}_{i=t+1}^{d-1}{(i-A_2)}^2{p}_i}{P_2}$

其中，P₁为C₁中像元数占水体指数C中总像元数的比例，P₂为C₂中像元数占水体指数C中总像元数的比例，i为像元灰度值，A为C的像元灰度平均值，A₁为C₁的像元灰度平均值，A₂为C₂的像元灰度平均值，p_i为灰度值为i的像元数占水体指数C中总像元数的比例。

2.根据权利要求1所述的基于多光谱影像的湖泊水储量遥感估算方法，其特征在于，所述步骤S1进一步包括：

收集和遴选遥感影像；

将收集和遴选的遥感影像进行辐射纠正和几何纠正；

根据纠正后的遥感影像提取反映湖泊水面分布特征的水体指数。

3.根据权利要求1所述的基于多光谱影像的湖泊水储量遥感估算方法，其特征在于，步骤S1中所述根据湖泊水体在多光谱遥感影像不同波段上的光谱响应特性提取反映湖泊水面分布特征的水体指数包括：

根据湖泊水体在Landsat MSS或HJ-1A/B影像不同波段上的光谱响应特性提取归一化差异水体指数；

或者，

根据湖泊水体在Landsat TM/ETM+影像不同波段上的光谱响应特性提取增强型归一化差异水体指数。

4.根据权利要求3所述的基于多光谱影像的湖泊水储量遥感估算方法，其特征在于，所述归一化差异水体指数NDWI采用下式计算：

$\mathrm{NDWI}=\frac{{\mathrm{\ρ}}_{\mathrm{Green}}-{\mathrm{\ρ}}_{\mathrm{NIR}}}{{\mathrm{\ρ}}_{\mathrm{Green}}+{\mathrm{\ρ}}_{\mathrm{NIR}}}$

其中，ρ_Green为Landsat MSS或HJ-1A/B影像上的绿光波段，ρ_NIR为Landsat MSS或HJ-1A/B影像上的近红外波段。

5.根据权利要求3所述的基于多光谱影像的湖泊水储量遥感估算方法，其特征在于，所述增强型归一化差异水体指数MNDWI采用下式计算：

$\mathrm{MNDWI}=\frac{{\mathrm{\ρ}}_{\mathrm{Green}}-{\mathrm{\ρ}}_{\mathrm{SWIR}}}{{\mathrm{\ρ}}_{\mathrm{Green}}+{\mathrm{\ρ}}_{\mathrm{SWIR}}}$

其中，ρ_Green为Landsat TM/ETM+影像上的绿光波段，ρ_SWIR为Landsat TM/ETM+影像上的短波红外波段。

6.根据权利要求1～5任意一项所述的基于多光谱影像的湖泊水储量遥感估算方法，其特征在于，步骤S3中所述根据湖泊等水位线数据集模拟湖泊水下地形包括：

根据湖泊等水位线数据集采用不规则三角网模拟湖泊水下地形，即湖泊水下地形由一系列互不重叠的不规则三角形组成的网状结构来表示，每个三角形的每个节点包括独立的经纬度坐标、水深值和坡度角；

步骤S4中所述根据模拟的湖泊水下地形和实测的水位数据计算不同时期的湖泊水储量包括：

根据模拟的水下地形和实测的水位数据将湖泊水体简化成一系列三棱柱的集合，整个湖泊水体的体积则为各三棱柱体积之和，计算公式为：

$V=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{S}_i\frac{h_i+h_{i+1}+h_{i+2}}{3}$

其中，n为三棱柱的个数，S_i为第i个三棱柱水体上表面的面积，h_i、h_i+1和h_i+2为三棱柱三条侧棱所在位置的水深值。

7.根据权利要求6所述的基于多光谱影像的湖泊水储量遥感估算方法，其特征在于，所述根据湖泊等水位数据集采用不规则三角网模拟湖泊水下地形包括：

根据湖泊等水位数据集采用ArcMap 9.3软件的三维分析工具中的TIN创建模块模拟湖泊水下地形。

8.根据权利要求6所述的基于多光谱影像的湖泊水储量遥感估算方法，其特征在于，所述湖泊水储量利用ArcMap 9.3软件三维分析工具中的面积-体积统计模块进行计算。