CN110335355A

CN110335355A - 一种大型浅水湖泊水面高自动计算方法

Info

Publication number: CN110335355A
Application number: CN201910638836.7A
Authority: CN
Inventors: 周航宇; 柯劲松; 杨智翔; 张冬; 张李荪; 刘金生
Original assignee: Jiangxi Water Conservancy Planning And Design Institute
Current assignee: China Railway Water Resources And Hydropower Planning And Design Group Co ltd; Jiangxi Water Conservancy Planning And Design Institute Co ltd
Priority date: 2019-07-16
Filing date: 2019-07-16
Publication date: 2019-10-15
Anticipated expiration: 2039-07-16
Also published as: CN110335355B

Abstract

本发明涉及测绘科学与技术，公开了一种大型浅水湖泊水面高自动计算方法，包括湖泊区域内DEM获取；湖泊区域内Landsat影像获取；坐标系统统一和投影；格网大小的选取；构建区域格网内水面面积‑水面高程关系曲线；获取区域格网内TM影像水域面积；区域格网内平均水面高的自动计算；自动粗差剔除和平差；相邻格网高程的比较；当相邻格网水面高程相差太大时，细分格网，重复格网大小的选取；当相邻格网水面高程相差符合要求时，则输出各格网水面高程。本发明，融合卫星瞬时遥感影像和实测DEM数据，在忽略浅水湖泊地形变化的情况下，通过自动格网划分和融合计算，可以获取格网内大湖瞬时水面高程，精度优于0.2米，解决了大型浅水湖泊水面瞬时高程获取难题。

Description

一种大型浅水湖泊水面高自动计算方法

技术领域

本发明涉及测绘科学与技术，具体是一种大型浅水湖泊水面高自动计算方法，特别是涉及遥感影像与实测DEM数据融合通过求取区域格网内的遥感水域面积和区域内水面面积-水面高程曲线自动计算水面高程的方法。

背景技术

在大型浅水型湖泊中，由于四周来水不定，且由于湖泊范围较大，导致上下游落差较大，在同一时刻，不同位置的进水口和出水口的水面高程均有较大差异，有的甚至相差十来米，造成湖水面呈极不规则状态，不同区域格网水面高程差异较大。

由于布设实时水位观测站费用较高，也就难以在大型湖泊布设足够多的实时水位观测站来满足水文计算需要，目前传统测量方法和观测方法难以获取某时刻全湖泊的瞬时水面高程，这也就难以定量获取湖泊容积，给水文计算和洪水分析带来不便。

发明内容

本发明的目的在于提供一种大型浅水湖泊水面高自动计算方法，以解决上述背景技术中提出的问题。

为实现上述目的，本发明提供如下技术方案：

一种大型浅水湖泊水面高自动计算方法，包括如下步骤：

（1）湖泊区域内DEM获取，先采用水下地形测量方法或者其他专业测量手段，获得湖泊内数字高程模型DEM；

（2）获取Landsat影像；

（3）采用ENVI专业图像处理软件，将Landsat影像图转换成与DEM数据同一椭球下相同投影方式下的影像图，以确保DEM和遥感影像在相同坐标系统和投影方式下；

（4）选取格网大小；

（5）构建区域格网内水面面积-水面高程关系曲线；

（6）获取区域格网内TM影像水域面积；

（7）自动计算区域格网内平均水面高；

（8）自动粗差剔除和平差；

（9）比较相邻格网高程；

（10）当相邻格网水面高程相差太大时，细分格网，重复（4）；当相邻格网水面高程相差符合要求时，则输出各格网水面高程。

作为本发明进一步的方案：步骤（5）中，区域格网内水面面积-水面高程关系曲线构建方法是：

（1）定义格网大小；

（2）根据DEM确定区域格网内水位Zs时对应的淹没面积；

（3）迭代过程或循环，根据不同的Zs，分别求取格网内低于或等于Zs高程的总面积；

（4）通过Zs及其对应面积，拟合构建格网内水面面积-水面高程关系曲线。

作为本发明进一步的方案：步骤（6）中，区域格网内TM影像水域面积计算方法是：

（1）对获取的遥感影像进行格网化处理，编写软件读取Landsat对应波段的灰度值，然后对像素灰度值进行归一化处理，获取区间-数量关系曲线；

（2）对MNDWI值进行区域细分，将归一化后MNDWI值域区域[-1，1]分成N等分，每区间长度为2/N，计算值域区间长度内像素点个数；

（3）对MNDWI值区间像素点个数进行统计分析，可发现MNDWI区域值具有两个近似正态分布区域，由于此两个波段是对水域敏感的波段，可认定两个正态分布分别为水域和非水域，确定两个近似正态分布区域，即水域和非水域；

（4）在两个近似正态分布区域之间，通过两极大值的最佳阈值确定法，找到区域中阈值C，通过对阈值C的判断便可获取区域内对应的水面面积。

作为本发明进一步的方案：步骤（7）中，区域格网内平均水面高的自动计算方法是：通过区域格网内构建淹没面积-湖水位的关系曲线和遥感影像相应区域内的水域面积，通过平滑样条曲线插值，便能获取相应的水位高程，在此，我们主要考虑平面坐标的确定和区域格网内平均水面高的获取。

其具体操作步骤为：（1）针对格网内所有判断为水域的平面格网坐标，取其平面格网平均值，当判断水域像素点个数为n时，平面格网中心坐标为（Xi，Yi），且平面坐标计算公式如式（1）。

（1）

（2）通过获取的水域面积和区域格网内构建淹没面积~湖水位曲线，求取水面高程。如果遥感影像相应区域获取水面面积为A，构建淹没面积~湖水位曲线时有相应点Ai<A< Ai+1时，水面高程Z可通过公式（2）获得：

（2）

作为本发明进一步的方案：步骤（8）中，自动粗差剔除和平差方法是：

（1）对所求的格网内水面高程，利用湖四周布设的实时水文观测站数据，删除明显粗差点，对于大于实时最高水位高程和小于实时最低水位高程点，根据实际情况进行判断删除；

（2）采用曲面拟合方法或者分区曲面拟合方法对所有格网水面高程点进行拟合，得出各点拟合残差；

（3）剔除拟合残差大于三倍中误差的点，再次进行曲面拟合或分区曲面拟合，直到所有点残差均在三倍中误差以内；

（4）求取拟合系数；

（5）对所有格网均采用拟合系数求取相应水位高程。

作为本发明进一步的方案：当分区太大时，会造成相邻格网内水面高程差距太大。比较格网内相邻的水面高程，步骤（10）中，当相邻格网水面高程差小于0.2米时，格网不再细分；当相邻格网水面高程差大于0.2米时，将格网细分，重定义格网。

作为本发明进一步的方案：读取Landsat对应波段的灰度值包括Band2和Band5的灰度值。

作为本发明再进一步的方案：删除明显粗差点是根据历史最高水位和历史最低水位进行判断。

与现有技术相比，本发明的有益效果是：解决了大型浅水型湖泊瞬时水面高程获取问题，为水文计算提供了准确基础数据，该方法可以适应于地势较平坦的湖泊水面高程计算；在采用改进的归一化水体指数法计算水域面积时，提出了一种基于两个极大值的阈值确定方法，可根据格网内像素的灰度值确定阈值进而自动判断是否水域；提出了一整套大型湖泊水面高计算方法。该方法解决了传统方法不能同时获取大面积水面高程的问题。

附图说明

图1为一种大型浅水湖泊水面高自动计算方法的流程示意图。

图2为一种大型浅水湖泊水面高自动计算方法中区域格网内DEM部分数据示意图。

图3为一种大型浅水湖泊水面高自动计算方法中格网内水面高程与水面面积关系示意图。

图4为一种大型浅水湖泊水面高自动计算方法中归一化指数阈值的确定方法。

图5为一种大型浅水湖泊水面高自动计算方法中MNDWI阈值确定示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1～5，本发明实施例中，一种大型浅水湖泊水面高自动计算方法，包括如下步骤：

（2）获取Landsat影像；

（4）选取格网大小；

（5）构建区域格网内水面面积-水面高程关系曲线；通过对区域格网大小的定义，根据DEM确定区域格网内不同水位Zs时对应的淹没面积As，循环迭代计算所有Zs时对应的淹没面积As，然后通过样条曲线进行内插，进而根据计算点（Zs，As）和内插点（Zs，As）对拟合构建格网内水面面积-水面高程关系曲线。

（6）获取区域格网内TM影像水域面积；编写软件读取对水体具有明显差异光谱特征的Band2和Band5的灰度值，对像素灰度值进行归一化处理，获取区间-数量关系曲线。然后对MNDWI值进行区域细分，将归一化后MNDWI值域区域[-1，1]分成N等分，每区间长度为2/N，计算值域区间长度内像素点个数。通过对MNDWI值区间像素点个数进行统计分析，可发现MNDWI区域值具有两个近似正态分布区域，由于此两个波段是对水域敏感的波段，可认定两个正态分布分别为水域和非水域，进而根据基于两个极大值的阈值判断方法，求取区域格网内水域面积大小。

（7）自动计算区域格网内平均水面高；针对格网内所有判断为水域的平面格网坐标，取其平面格网平均值，当判断水域像素点个数为n时，平面格网中心坐标为（Xi，Yi），且平面坐标计算公式如式（1）。

（1）

通过获取的水域面积和区域格网内构建淹没面积~湖水位曲线，求取水面高程。如果遥感影像相应区域获取水面面积为A，构建淹没面积~湖水位曲线时有相应点Ai<A<Ai+1 时，水面高程Z可通过公式（2）获得：

（2）

（8）自动粗差剔除和平差；通过历史最高水位和历史最低水位对区域格网内水面高程进行判断，删除粗差点，然后采用分区曲面拟合方法对所有点进行拟合计算，删除拟合残差大于三倍中误差的点，然后循环拟合，直到所有格网点残差均在三倍中误差之内，此时求解拟合系数，并计算相应水面高程，进而输出水面高程点坐标。

（9）比较相邻格网高程；

步骤（1）:

定义鄱阳湖范围，通过专业测量手段获取湖泊内数字高程模型，可采用单波束测深仪+GPS测量水下地形，航摄与Lidar摄影测量方法测量陆地地形，根据相关软件生成数字高程模型。

步骤（2）:

根据（1）定义的范围，下载相应的Landsat影像，每景影像必须覆盖鄱阳湖区域；

步骤（3）:

Landsat影像采用UTM投影，DEM坐标是2000坐标系，采用高斯-克吕格投影。由于两者中央经线上长度比不一样，需经过专业软件进行转换。本发明采用ENVI专业图像处理软件，将Landsat影像图转换成与DEM数据同一椭球下相同投影方式下的影像图（坐标与DEM数据一致），以确保DEM和遥感影像在相同坐标系统和投影方式下。

步骤（4）:

格网大小的选取，为了能够自动确定格网大小，逐渐满足需要。格网大小改变需要按照一定的规律进行。首先将区域格网大小确定为一幅1:1万地形图，当不能满足要求时，将图幅分为k×k个格网， k值分别取为1,2,3…

步骤（5）:

① 格网的定义

格网化是研究地理空间信息分布的一种重要方法。在涉及地理特征的人口分布、建筑物空间分布、重力异常分布、GDP分布等很多领域，经常采用格网化研究方法。目前采用的地形图国际分幅是指按国际统一地形图分幅原则划分的地理图幅，是经纬网梯形分幅法，即按规定的经差和纬差划分图幅，以经纬度格网框架为基础进行的，也是格网化的一种具体应用。例如：国际基本比例尺地形图分幅与编号，以1:100万地形图为基础，按规定的经差和纬差划分图幅。从赤道起算，每纬差4°为一行，至88°，南北半球各分为22横列，依次编号A、B、... V；由经度180°西向东每6°一列，全球60列，以1-60表示，如海南所在1：100万图在第5行，第49列，其编号为 E49 ；北京为J50，例如：每幅1：100万地形图划分为96行96列，按经差3′45″纬差2′30″分成9216幅1：1万地形图。在世界各国的小比例尺地形图分幅中，大多以经纬度格网框架为基础，生产出来的地形图、影像图、数字高程模型也按照经纬度来划分存储。

在库（湖）容计算中，由于面积大，所涉及到的地图数据多，一次性读入，不仅造成计算机运行时间长、速度慢，有时甚至出现难以运行、发生死机现象。为提高效率，减少一次性读入数据量，本文在水位推算和湖容计算过程中，均考虑仅读入单幅数据，并可根据需要选择对单幅地图数据再进行格网化。划分格网的原则是当相邻区域格网推导计算的水位高小于某一规定值时，此时的格网无需继续细分。

② 根据DEM确定区域格网内水位Zs时对应的淹没面积

数字高程模型（DEM）由等间隔海拔数据排列组成。区域格网是以单幅小比例尺地图为基础进行格网划分的，可知，区域格网一定是梯形，而大多数DEM格网都是正方型（X和Y方向的分辨率相同）。如附图2。

在采用DEM计算区域格网内水位Zs对应的淹没面积时，应考虑以下几点：

① DEM格网点中心坐标是否落在湖区范围线内；

② DEM格网点中心坐标是否落在区域格网内；

③ 区域格网不宜过小，应远大于DEM格网，否则需考虑范围线和区域格网外围线与DEM格网相交处的实际格网面积；

④ 当区域格网远大于DEM格网时，区域内所有低于Zs的高程格网点个数对应的面积可采用公式（3）计算：

As=n*D<sup>2 </sup>

（3）

N代表区域格网内低于Zs高程所有像素数目，D表示DEM格网大小（通常每一个格网大小是一个平方数，例如x和y方向的分辨率是相同）

③ 迭代过程或循环。

对区域格网内所有水位Zs高程依次求取对应面积，Zs确定方法有两种：

Zs可以根据精度要求进行确定，比如：当DEM精度为0.1米时，可以从区域格网所有DEM中先读出最小高程和最大高程，最小高程与最大高程中间逐渐累加0.1米与区域所有DEM格网高程进行比较，利用（3）式计算对应面积。如果增加若干米后，计算对应面积与不增加前相等，此数据可忽略。例如：当小于等于22.3米和小于等于22.5米时的对应淹没面积相等时，前面才是正确淹没面积。

读取所有DEM格网高程数据，然后对高程数据进行排序，按照排序后的高程依次确定对应淹没面积。上述格网分辨率为10米，可知在有效区域范围内的格网个数和面积如表1。

表1 DEM高程-面积对应表

④ 拟合构建格网内水面面积~水面高程关系曲线

通过第三步可获取格网内水面面积~水面高程对应关系。以水面高程为纵坐标，水面面积为横坐标，可以构建格网内水面面积~水面高程关系曲线。可采用解析表达式逼近离散数据的方法和最小二乘法来进行拟合。一般采用最小二乘原理对曲线进行拟合，常用曲线有：直线、二次多项式、三次多项式、三次样条插值等。当采用曲线拟合的方法进行计算时，由于格网内水面面积~水面高程关系点较多，而当采用某个函数进行逼近时，会将曲线点上所有点进行平滑拟合，加剧了局部点拟合的误差，所以不宜采用某个函数来模拟定义区间的所有样本，宜采用分段函数进行逐一模拟。当DEM精度较高时，可采用简单的样条曲线进行内插，本研究采用样条曲线进行内插，如图3。

步骤（6）:

对获取的遥感影像进行格网化处理，利用软件读取Landsat对应波段的灰度值，此处读取Band2和Band5的灰度值，对像素灰度值进行归一化处理，获取区间-数量关系曲线。然后对MNDWI值进行区域细分，将归一化后MNDWI值域区域[-1，1]分成N等分，每区间长度为2/N，计算值域区间长度内像素点个数。通过对MNDWI值区间像素点个数进行统计分析，可发现MNDWI区域值具有两个近似正态分布区域，由于此两个波段是对水域敏感的波段，可认定两个正态分布分别为水域和非水域，在两个近似正态分布区域之间，利用下列方法，可找到区域中阈值C，通过对阈值的判断便可获取区域内对应的水面面积，确定阈值C的方法如图4。

1.将MNDWI值域区域归一化到区域[-1，1]；

b.由各种研究数据可知，经归一化后的阈值一般处于[-0.2，0.2]之间，为了更好的确定阈值，将[-1,1]区间分为2000等分或者更多，设置Number Range=2000。

c.利用改进的归一化水体指数法计算每个像素的MNDWI因子，将计算的因子按区间统计数量，可得到区间～数量统计图。利用程序，可考虑值域的相似性，例如：如果区间分为10000，但每个区间上数量均不大于100，此时可认为区间分得过细，宜设置值域不同数量数大于区间数量的一半，通过设置循环程序可完成此区间个数值的确定。

d.为简化处理，我们将区间 [k,k+1]所对应值域数量Y设定为在k时所对应的数量Y，即得到（-0.3+0.6/60*k, Y）点序列；

e.由以上图区间数量统计图可知，阈值的确定转化为寻找两个极大值区间内极小值的过程，如图5；

f.利用第一个极大值与另一个极大值相隔区间数至少大于等于4，可以完成两个极大值的查找。

g.建立与上述点序列一致的临时点序列Temp（k,-0.3+0.6/60*k, Y）；

h.对所有值域进行循环，找到函数值域最大值Yk1，此时得到对应的K1点（MNDWIk,Yk），令MaxNumber0=k1,MaxNumber1=k1；

i.在临时点序列中删除MaxNumber1最大值，并查找次大值，即查找删除最大值后的点序列中Y最大值K2；

j.如果abs(K1-K2)的区间数小于等于4，则说明K2仍在第一个最大值所属区间内，此时，令MaxNumber1=k2，重复第7步；

如果abs(K1-K2)的区间数大于3，则找到第二个峰值MNDWI_k2，此时将临时点序列中的MNDWI_k2与原序列相比；当MNDWI_k2与原序列中K3相一致时，此时得到原序列第二个峰值序号MaxNumber2；

k.第一个峰值序号是MaxNumber0；

l.通过查找两个极大值之间的极小值，可以确定相对应的阈值；

m.通过判断阈值获取区域格网内水体面积。

步骤（7）:

（1）针对格网内所有判断为水域的平面格网坐标，取其平面格网平均值，当判断水域像素点个数为n时，平面格网中心坐标为（Xi，Yi），且平面坐标计算公式如式（1）。

（2）通过获取的水域面积和区域格网内构建淹没面积~湖水位曲线，求取水面高程。如果遥感影像相应区域获取水面面积为A，构建淹没面积~湖水位曲线时有相应点A_i<A<A_i+1时，水面高程Z可通过公式（1）获得。

步骤（8）:

（1）对所求的格网内水面高程，利用湖四周布设的实时水文观测站数据，删除明显粗差点（根据历史最高水位和历史最低水位进行判断），对于大于实时最高水位高程和小于实时最低水位高程点，根据实际情况进行判断删除。

（4）求取拟合系数；

（5）对所有格网均采用拟合系数求取相应水位高程；

步骤（9）:

循环调用相邻格网高程，并进行比较；

步骤（10）:

当相邻格网水面高程相差太大时，细分格网，重复（4）；当相邻格网水面高程相差符合要求时，则输出各格网水面高程。

本发明的工作原理是：大湖水面由于存在进水口、出水口水位高程的差异，在较远距离时，水面高程也必然存在差异。当距离缩小到一定区域范围时，由于水位平衡的快速调整（如果忽略风和波浪带来的变化），此时，格网区域内的水位可以看成一个静止的平面，当格网区域内地势平坦，且水域面积与陆地面积共存时，区域内一定的水域面积对应区域内一定的水位高程；同样，一定的水位高程必定对应一定的水域面积。根据水面具有延续性的特点，本发明在自动调整区域格网大小的基础上，通过融合区域格网内瞬时遥感影像和实测DEM，根据区域格网内水域面积大小和水面面积-水位高程的关系曲线，求取该格网瞬时水面高程。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims

1.一种大型浅水湖泊水面高自动计算方法，其特征在于，包括如下步骤：

（2）获取Landsat影像；

（4）选取格网大小；

（5）构建区域格网内水面面积-水面高程关系曲线；

（6）获取区域格网内TM影像水域面积；

（7）自动计算区域格网内平均水面高；

（8）自动粗差剔除和平差；

（9）比较相邻格网高程；

2.根据权利要求1所述的一种大型浅水湖泊水面高自动计算方法，其特征在于，步骤（5）中，区域格网内水面面积-水面高程关系曲线构建方法是：

（1）定义格网大小；

（2）根据DEM确定区域格网内水位Zs时对应的淹没面积；

3.根据权利要求1所述的一种大型浅水湖泊水面高自动计算方法，其特征在于，步骤（6）中，区域格网内TM影像水域面积计算方法是：

（3）对MNDWI值区间像素点个数进行统计分析，确定两个近似正态分布区域，即水域和非水域；

4.根据权利要求1所述的一种大型浅水湖泊水面高自动计算方法，其特征在于，步骤（7）中，区域格网内平均水面高的自动计算方法是：通过区域格网内构建淹没面积-湖水位的关系曲线和遥感影像相应区域内的水域面积，通过平滑样条曲线插值，便能获取相应的水位高程。

5.根据权利要求1所述的一种大型浅水湖泊水面高自动计算方法，其特征在于，步骤（8）中，自动粗差剔除和平差方法是：

（4）求取拟合系数；

（5）对所有格网均采用拟合系数求取相应水位高程。

6.根据权利要求1所述的一种大型浅水湖泊水面高自动计算方法，其特征在于，步骤（10）中，当相邻格网水面高程差小于0.2米时，格网不再细分；当相邻格网水面高程差大于0.2米时，将格网细分，重定义格网。

7.根据权利要求3所述的一种大型浅水湖泊水面高自动计算方法，其特征在于，读取Landsat对应波段的灰度值包括Band2和Band5的灰度值。

8.根据权利要求5所述的一种大型浅水湖泊水面高自动计算方法，其特征在于，删除明显粗差点是根据历史最高水位和历史最低水位进行判断。