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 栅格地图数据矩阵 （上图excel设计底稿下载)

当要求堰塞湖水位海拔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 平台的优势就是使用它的数据存储基础服务，用户关心的空间范围和显示服务，把所有精力都集中在算法本身上面。对地址空间的处理操作并不提供现成的服务。

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

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