当前位置: 代码迷 >> GIS >> Grass Gis 崖谷堰塞湖预测模型r.lake剖析-湖心灌水洪泛算法
  详细解决方案

Grass Gis 崖谷堰塞湖预测模型r.lake剖析-湖心灌水洪泛算法

热度:393   发布时间:2016-05-05 06:05:27.0
Grass Gis 峡谷堰塞湖预测模型r.lake剖析--湖心灌水洪泛算法

Grass Gis 是20世纪80年代初,美国军方建筑工程研究实验室(USA/CERL)的Bill Gorgan负责遴选一款具有土地管理、环境规划、环境评估的GIS软件。它有400多个地质分析模块,今天我们剖析一下r.lake算法。


1. 堰塞湖预测建模r.lake能做什么

如下图1左为一张地形栅格图,也是最常见的2D平面图。从图中可以看到许多的凹地盆地,当有山洪暴雨时就会形成堰塞湖。 Grass Gis 的r.lake就是能够预测堰塞湖形状大小的模块。例如,在图1左所示的红色圆圈位置开始注水当水位海拔高度为113.4米时,形成的堰塞湖的位置形状如图1右所示。Grass可以形成动画动态展示不断的提高水位后堰塞湖形态的变化。


我国长江三峡大坝水位设计海拔185米,蓄水库区淹没129座城镇,淹没耕地1.94万公顷东起宜昌,西至重庆662.9Km沿长江两岸分水岭范围。在设计时就需要对185米高湖所能淹没的地区有非常精确的预测,如果沿岸有缺口造成计划外的地区被淹没就需要处理缺口。


Grass Gis 的r.lake 就可以完成这个事情,它的核心算法就是“湖心灌水洪泛算法“,它模拟了水在一点不断涌出时淹没了周边地区。



图1 堰塞湖预测模型演示


2. 数据降维简化模型

当一个点开始蓄水时,周边山岭的高度和形态都可以阻挡洪水的进一步蔓延,就像一个围栏一样。那怎样才能根据水位的海拔高度找出这些围栏的位置?首先想到的是三维建模以及空间关系计算等复杂的方式,其实不然。所有复杂的三维空间计算都要降维到二维平面计算,简化计算。


一个地区的遥感影像被处理成高程图也就是栅格图后,就变成了一个标准的M*N的数字矩阵。矩阵的每个点对于地图上的一个矩形区域,该点的值就是该区域的平均高度。图1就是这样的一个二维栅格图,只是每个像素被赋予不同的灰度颜色展现出来。一个复杂的三维地势图化简为矩阵后就变成了如图2a所示的二维数组,图中的数字为某块土地的海拔。


图2 栅格地图数据矩阵


当要求堰塞湖水位海拔9米时,满足条件的区域如图2b所示的三块区域。这样问题好像看起来非常简单了, 就是把数组遍历一遍找出<=9的区域不就行了?但是复杂的问题来了,当洪水爆发点在 6C 点,那9米的水只能在红色区域聚集成堰塞湖,黄色和蓝色都不应该被淹没,因为水被红色区域周围20多米高的地势挡住了。


现在的问题已经被简化成了一个简单的数据结构算法了。以 6C 点逐渐向往扩张,就像湖心有洪水不断涌出直到把周围的相邻8个格子中低于9米的区域淹没。6C为湖心的周边5B,5C,5D,6B,6D,7B,7C,7D被淹没了。接着以它们为湖心重复相同的操作,尝试淹没周围的8个邻居区域。不断重复直到没有新的可淹没区域被找到。这个算法很简单,就不给出代码了。


3. 以一个或多个堰塞湖区为洪水点

 Grass Gis 的r.lake提供了更为复杂的建模。不是以某个点为洪水发生地,而是以某个湖区在涌出洪水,甚至多个湖区一起涌出。这样每次都以上一轮模拟输出的湖区为下次模拟计算的输入,不断增高湖水高度,达到动画效果输出。


这个算法思想和上面单点的算法相似,只是每轮需要扫描M*N个像素点,每次是旁边8个区域只要有一个有水的点就试图向自己淹没,直到一轮没有新的扩张点就停下来。

</pre><p><pre name="code" class="cpp">int pases = (int)(rows * cols) / 2;int  water_level = .... // 输入的要求的湖面海拔高度int *in_terran; // row * col 输入的栅格地图,数值为地势的海拔高度int *in_wanter; // row * col 输入的已经有水的湖区矩阵. 0无水,>0 为湖水平面到湖底的深度int *out_water; // row * col 输出计算完毕的湖区矩阵. 初始为0无水,>0 为湖水平面到湖底的深度copy_seed_water(in_wanter, out_water); // 初始输出湖区为输入湖区int water_window[3][3]; // 某点为中心的3*3矩阵窗口,该窗口沿着地图每个点移动一遍为一趟for (pass = 0; pass < pases; pass++) {	curcount = 0;	/* 每个点扫描 */	for (row = 0; row < rows; row++) 	{	    for (col = 0; col < cols; col++) 		{			//该点为中心去3*3矩阵			load_neibour_window_values(out_water, water_window, rows, cols, row, col); 			//该点周围8个邻居是否有水往自己这边淹没			if (is_near_water(water_window) == 1) 			{				if (in_terran[row][col] < water_level) 				{					// 该地势低,被淹没					out_water[row][col] =  water_level - in_terran[row][col];					curcount++;				 }				 else 				 {					out_water[row][col] = 0;	/* 该地势比湖水面高 */				 }			}	    }	}	if (curcount == lastcount)	    break;		/* 没有发现新的淹没点,完成. */	lastcount = curcount;	curcount = 0;	}}return out_water;


4. Grass Gis平台提供的服务支持

上面的算法没有看到任何Grass Gis 平台API的使用。其实栅格数据的输入,输出和图形化显示都是透明的交给了grass 平台了。向湖区的梯度渐变颜色表的建立就是调用API并绑定湖区地图的,但是开发人员并不知道该颜色表是怎样存储并绑定的。它提供的功能包括:

  • 各种地质地图数据文件在内存中的直接访问
  • 处理完成的地质数据向grass数据库的内存直接写
  • 地质地图颜色表的生成并绑定
  • 多个地质地图的三维叠加显示
  • 统一管理用户设定的的数据空间范围
所以使用Grass Gis 平台的优势就是使用它的数据存储基础服务,用户关心的空间范围和显示服务,把所有精力都集中在算法本身上面。对地址空间的处理操作并不提供现成的服务。

版权声明:本文为博主原创文章,未经博主允许不得转载。

版权声明:本文为博主原创文章,未经博主允许不得转载。



版权声明:本文为博主原创文章,未经博主允许不得转载。