元胞自动机-森林火灾模拟

引入：

    元胞自动机，英文名及缩写：cellular automata，CA。最初是由冯诺依曼在二十世纪五十年代为模拟生物自保的自我复制而提出的，但是当时并未受到重视。后来才逐渐发展起来，著名的“生命游戏”就是其一个应用案例。但今天想和大家分享的案例并不是生命游戏，而是来自近期课程实验中的内容：根据元胞自动机原理编程实现“森林火灾模拟”。

---

正文:

首先还是从元胞自动机的基本理论说起吧，先谈一下什么是元胞自动机。

元胞自动机的定义

    ①物理学定义： 元胞自动机是定义在一个由离散、有限状态的元胞 组成的元胞空间上，按照一定的局部规则，在离散时 间维度上演化的动力学系统。
    ①数学定义：设d表示空间维数，k表示元胞的状态，并在一个有限集合S中取值，r表示元胞的邻居半径。Z是整数集，表示一维空间，t表示时间。元胞自动机的动态演化就是在时间上状态组合的变化，可记为:

F:S(t)Z->S(t+1)Z

可以看到，上述定义中提到了以下几个关键词：*元胞、元胞空间、邻域、规则*。下面结合一张图片来堆对其进行简单解释。

(1)元胞：是元胞自动机的基本单元，每一个元胞都有记忆存储的功能，并且元胞自动机系统中的所有元宝都按照动力规则不断更新。简单来说，元胞就类似上述棋盘图片中的棋子，而每个元胞的记忆则用于描述与棋子自身的位置与颜色。
(2)元胞空间：是指元宝分布的空间网格的集合，分为一维、二维、三维元胞空间，比如上述图片中的棋盘空间就是一个二维元胞空间。当然，此次模拟的森林火灾对应的也是一个二维元胞空间。
元胞和元胞空间用于表示系统的静态成分，就是说：仅仅只有这两者，是无法支撑元胞自动机的正常运行的，为此，我们需要引入系统的动态组分，即下面要介绍的演化规则。
(3)邻域：元胞自动机的演化规则是定义在局部范围内的，即：一个元胞在（t+1）时刻的状态决定于该元胞在 t 时刻本身的状态和它的邻居元胞的状态。
二维空间中的邻域常见的有：VonNeumann邻居、Moore邻居、扩展Moore邻居等。
此次要用到的就是VonNeumann邻域。

(4)规则：元胞自动机根据规则进行局部元胞间的相互作用，而引起全局变化，规则支配着元胞自动机的整个动力学行为。即：根据元胞当前状态及其邻居状态确定下一时刻该元胞状态。

---

森林火灾模拟问题规则定义与分析

    好了，有了上面的基本概念，我们先来做一下分析。上面我们提到，此次森林火灾模拟是针对二维元胞空间，为元胞的动态演变定义局部演变规则，最终实现元胞的自动演变。

森林火灾模拟问题简析：

可将森林视为元胞自动机模型中的二维元胞空间，而森林中的树木视为元胞空间中(x,y)坐标处的一个元胞，具备一定的状态；下面来看一下该系统中元胞的运动规则：

---

※设计森林火灾的模拟的规则—— 森林灾模拟时元胞有3个不同的状态（分别为状态0：空位；状态2：燃烧着的树木；状态1：树木），采用的邻域为VonNeumann邻居。 演化规则如下：
（1）若某树木元胞的4个邻居中有燃烧着的，那么该元胞下一时刻的状态是燃烧；
（2）一个燃烧着的元胞，在下一时刻变成空位；
（3）所有树木元胞（状态为1）以一个低概率开始燃烧（模拟闪电造成森林火灾）；
（4）所有空位元胞以一个低概率变为树木（模拟新树木的生长）。

---

森林火灾模拟邻域及规则阐述

①邻域：规定使用VonNeumann邻居，那么一个位置在(x,y)处的元胞的邻域范围Δ中其他元胞坐标为：
【1】(x-1,y);【2】(x,y-1);【3】(x,y+1);【4】(x+1,y) ②规则：对应一个元宝的邻域范围Δ，若：
【1】第一条规则：搜寻邻域空间Δ的四个元胞状态，若有1个为燃烧状态'2'，则：(x,y)位置的元宝状态在下一刻状态也变为：燃烧'2';
【2】第二条规则：若(x,y)中心元胞位置处状态为：燃烧状态，则下一刻变为：空地’0‘；<>br 【3】第三条规则：给定一个低概率使状态为1的未燃烧树木元胞下一刻开始燃烧，状态2；
【4】第四条规则：若(x,y)中心元胞位置处状态为：空地，则下一刻在一定的给定概率条件下变为：树木’1‘；

---

森林火灾模拟问题求解步骤及其图形库：

---

求解步骤

基于上述分析，将该问题的求解步骤描述为：
①初始化格网矩阵，作为元胞空间来存放元胞群体；
②初始化占据网格总数一定比例的元宝群体，用于表述状态为1的正常树木；
③绘制元胞群体到图形界面上，用绿色标识状态为1的正常树木；用黑色表示状态为0的空地；用红色标识状态为2的燃烧中的树木；
④根据规则获取元胞的下一时刻状态；
⑤计算新生树木中发生火灾的个体数目。预先给出新生状态为1的树木发生火灾的概率，在此前提条件下计算并获取随机发生火灾后的元胞状态；
⑥重新绘制元胞群体到屏幕上。
⑦循环执行④-⑥，中间可以使其休眠1.5s，以便于观察结果的变化过程。
每个过程都被定义成一个函数，见后面代码部分有详细注释。

图形库介绍

编程过程中用到了EasyX图形库。
EasyX 是针对 C++ 的一个轻量级图形库，可以帮助 C 语言初学者快速上手图形和游戏编程。提供了绘图函数、鼠标与键盘事件响应函数、以及资源文件读取函数等。同时针对绘图部分也提供了批量绘图函数，解决了窗体刷新过程中出现的窗口闪烁的问题，这在下面的代码部分都有体现。

---

森林火灾模拟问题-代码[C语言]：

先贴个运行结果截图

#include<stdio.h>
#include<stdlib.h>
#include<graphics.h>
#include <time.h>
#include <math.h>

//定义树木状态->四邻域
#define no_tree 0	//空位
#define burning 2	//燃烧
#define nor_tree 1	//树木

//定义窗体尺寸
#define WIDTH 600	//宽度
#define HEIGHT 600	//高度
#define treeSize 4	//树木尺寸
#define treeNum	150	//树木数量
#define Row		150	//行数
#define Col		150	//列数
#define burning_rate	 0.0001//初始燃烧概率:0.01%
#define proTree_rate     0.02//空地长出新树木的概率：10％
#define unsafe_rate		 0.00005//在生长出新树木期间再次失火的概率：0.005%
#define N 10000

//定义位置结构体
struct POS
{
	int x;
	int y;
};

void InitBoard();//初始化布局网格
void InitTree(int treeState_Matrix[treeNum][treeNum]);//初始化树木种群
void drawTree(int treeState_Matrix[treeNum][treeNum]);//绘制树木种群
void InitRandTreePos();//生成单个火源初始位置
void getNextTreeState(int treeState_Matrix[treeNum][treeNum]);//获取森林树木的下一时刻状态
int GetArea_Tree(int i, int j); //检测VonNeumann邻域内是否存在燃烧树木
int isRight(float state);//产生一个随机数，用于判断树木下一时刻是否会燃烧/长出新树木
void burn_again();//在长出新树木时再次失火

/*
演化规则如下：
	（1）若某树木元胞的4个邻居中有燃烧着的，那么该元胞下一时刻的状态是燃烧；
	（2）一个燃烧着的元胞，在下一时刻变成空位；
	（3）所有树木元胞（状态为2）以一个低概率开始燃烧（模拟闪电造成森林火灾）；
	（4）所有空位元胞以一个低概率变为树木（模拟新树木的生长）。
*/

int treeState_Matrix[treeNum][treeNum];//树木状态矩阵

//主函数：程序入口
int main(int args,char* argv)
{
	initgraph(WIDTH, HEIGHT,SHOWCONSOLE);//初始化窗体：500*500
	HWND hwnd = GetHWnd();//获取窗体句柄
	SetWindowText(hwnd,_T("林森火灾模拟"));//设置窗体标题
	SetWindowPos(hwnd, NULL, 600, 200, HEIGHT, WIDTH, 0);//设置窗体位置
	//1、初始化网格布局
	InitBoard();
	srand((int)time(NULL));//随机数种子
	//2、初始化数目状态矩阵
	InitTree(treeState_Matrix);
	//3、绘制树木
	drawTree(treeState_Matrix);
	//4、初始化燃烧树木的位置
	InitRandTreePos();
	//5、绘制森林-树木
	drawTree(treeState_Matrix);
	
	//开启批量绘图模式
	BeginBatchDraw();
	while (true)//进入循环迭代
	{
		//6、获取下一时刻的树木状态矩阵
		getNextTreeState(treeState_Matrix);
		//7、仍然有可能产生新的火源
		burn_again();
		//8、重新绘制森林-树木
		drawTree(treeState_Matrix);
		//显示绘图结果
		FlushBatchDraw(); // 显示绘制 
		//停顿0.5s
		Sleep(500);
	}
	EndBatchDraw();//结束批量绘图模式
	closegraph();//关闭设备窗口
	return 0;
}

/*1、初始化布局网格*/
void InitBoard()
{
	setlinecolor(RGB(0, 0, 0));//黑色线条
	for (int i = 0; i < treeNum; i++)
	{
		line(0, i*treeSize, WIDTH, i*treeSize);
		line(i*treeSize, 0, i*treeSize, HEIGHT);
	}
}


/*2、初始化树木种群*/
void InitTree(int treeState_Matrix[treeNum][treeNum])
{
	int i, j,count=0;
	float randVal;
	for ( i = 0; i < Row;i++)
	{
		for (j = 0; j < Col;j++)
		{
			//随机概率大于0.4，生成树木：状态2
			randVal = rand() % (N + 1) / (float)(N + 1);
			if (randVal>=0.40)
			{
				 treeState_Matrix[i][j] = 1;//[i,j]表示ROW*COL规格的网格处对应的一个树木：产生树木
				 count++;
 			}
 			else
  			{	
 				treeState_Matrix[i][j] = 0;//否则为空位
 			}//end if
		}//end for内循环
	}//end for外循环
	printf("初始化正常树木状态矩阵完毕-%d\n",count);
}

/*3、绘制树木种群*/
void drawTree(int treeState_Matrix[treeNum][treeNum])
{
	int i, j;
	for (i = 0; i < Row; i++)
	{
		for (j = 0; j < Col; j++)
		{
			if (treeState_Matrix[i][j]==1)//正常树木：状态1
			{
				setfillcolor(RGB(0, 255, 0));//设置填充颜色：绿色
				fillrectangle(i*treeSize, j*treeSize,
					i*treeSize + treeSize, j*treeSize + treeSize);
			}else if (treeState_Matrix[i][j]==2)//燃烧的树木：状态2
			{
				setfillcolor(RGB(255,0,0));//设置填充颜色：红色
				fillrectangle(i*treeSize, j*treeSize,
					i*treeSize + treeSize, j*treeSize + treeSize);
			}
			else if (treeState_Matrix[i][j]==0)//空地：状态0
			{
				setfillcolor(RGB(0,0,0));//设置填充颜色：黑色
				fillrectangle(i*treeSize, j*treeSize,
					i*treeSize + treeSize, j*treeSize + treeSize);
			}
		}
	}
}

/*4、生成初始燃烧位置,概率：buring_rate*/
void InitRandTreePos()
{
	int i, j, count = 0;
	float randVal;
	for (i = 0; i < Row; i++)
	{
		for (j = 0; j < Col; j++)
		{
			if (treeState_Matrix[i][j]==1)
			{
				//随机概率大于0.4，生成树木：状态2
				randVal = rand() % (N + 1) / (float)(N + 1);
				if (randVal < burning_rate)
				{
					treeState_Matrix[i][j] = 2;//[i,j]表示ROW*COL规格的网格处对应的一个树木：产生树木
					count++;
				}//end if
			}
		}//end for内循环
	}//end for外循环
	printf("初始化燃烧树木状态矩阵完毕-%d\n",count);
}

/*5、获取森林树木的下一时刻状态*/
void getNextTreeState(int treeState_Matrix[treeNum][treeNum])
{
	int i, j,state=0;
	for (i = 1; i < Row-1; i++)
	{
		for (j = 1; j < Col-1; j++)
		{
			state = treeState_Matrix[i][j];//记录当前树木状态
			if (state==1)//如果为正常树木:1
			{
				if (GetArea_Tree(i, j))//判断下一刻是否会被引燃，如果返回1[表示VonNeumann邻域内存在燃烧树木]
				{
					treeState_Matrix[i][j] = 2;
				}
			}//end if-1
			else if (state==2)//如果为燃烧树木:2
			{
				treeState_Matrix[i][j] = 0;//下一刻变为空地
			}//end if-2
			else if (state==0)//如果为空地
			{
				if (isRight(proTree_rate))//判断是否有新树木的产生
				{
					treeState_Matrix[i][j] = 1;
				}
			}//end if-3
		}//end if-内
	}//end for-外
	printf("下一时刻状态获取完毕\n");
}

/*6、在长出新树木时有一定概率再次失火*/
void burn_again()
{
	int i,j,state = 0;
	float randVal = rand() % (N + 1) / (float)(N + 1);
	if (randVal <= unsafe_rate)//在失火概率范围内
	{
		for (i = 1; i < Row - 1; i++)
		{
			for (j = 1; j < Col - 1; j++)
			{
				state = treeState_Matrix[i][j];//记录当前树木状态
				if (state==1)//正常树木
				{
					if (isRight(burning_rate))
					{
						treeState_Matrix[i][j] = 2;
					}
				}
			}
		}
	}
}

//检测VonNeumann邻域内是否存在燃烧树木
int GetArea_Tree(int i, int j)
{
	int sum = 0;
	if ((treeState_Matrix[i - 1][j] == 2) ||
		(treeState_Matrix[i][j - 1] == 2) ||
		(treeState_Matrix[i][j + 1] == 2) ||
		(treeState_Matrix[i + 1][j] == 2)
		)
	{
		return 1;
	}
	return 0;
}

//产生一个随机数，用于判断树木下一时刻是否会燃烧
int isRight(float state)
{
	float randVal = rand() % (N + 1) / (float)(N + 1);
	if (randVal <= proTree_rate)//在生长树木概率范围之内
	{
		return 1;
	}
	return 0;//不//在燃烧概率范围之内
}

结尾:

    关于元胞自动机，再小补充一点copy过来的感受。元胞自动机不同于一般的动力学模型，因为元胞自动机不是由严格定义的物理方程或函数确定的，而是用一系列模型构造的规则构成。凡是满足这些规则的模型都可以算作是元胞自动机模型。
    因此，元胞自动机是更像一类模型的总称，或者说是一个方法框架。其特点是时间、空间、状态都离散，每个变量只取有限多个状态，且其状态改变的规则在时间和空间上都是局部的。而这些变量在局部时空上的改变，也正是由系统运行的核心规则所驱动的，在不同的时刻，展现出特定的表现型。穿插在这个森林火灾模拟案例中的这种编程思想，类似于时态GIS数据模型中的[基态修正模型](https://malagis.com/spatial-databases-103-temporal-gis-data-model.html)[感兴趣的朋友可以点击超链接看一下]。
    因为刚接触元胞自动机的基础知识，对其了解并不多，也只能套用之前学过的概念来对其进行理解，如果有出错的地方，请大家多多指教。另外，有关生命游戏和初等元胞自动机的内容有机会的话下次再和大家分享，谢谢！