Моделирование любого природного процесса — всегда довольно сложная задача, особенно если речь идет о погодных и климатических явлениях. По сути, так называемая процессно-ориентированная модель всегда опирается на некоторый набор допущений, позволяющих относительно просто описать определенную среду. Эти допущения являются необходимостью, поскольку мы не можем учесть каждый процесс или траекторию движения каждого объекта/молекулы/частицы. Поэтому «под капотом» таких моделей всегда находятся знания о реальном мире, обернутые в набор уравнений/правил, с некоторой точностью аппроксимирующих окружающую среду.
Я не являюсь экспертом в области моделирования на основе процессов, однако недавно я наткнулся на увлекательную статью Кристиана Хилла. Считаю, что это очень добротное и доступное введение в моделирование в условиях реального мира, которое может понять любой человек, обладающий базовыми математическими знаниями. Сегодня я попытаюсь изложить свою точку зрения на концепцию, описанную в вышеупомянутой статье, и построить чуть более сложную и реалистичную модель.
Код этой статьи вы можете найти на GitHub. Код, который вы увидите ниже, полностью написан мной, я лишь воспользовался подходом Кристиана.
Клеточный автомат
Для построения симуляции лесного пожара мы будем опираться на концепцию клеточного автомата (КА), который представляет собой вычислительную модель, используемую для моделирования сложных систем и процессов. КА состоит из N-мерной сетки ячеек.
- Каждая ячейка имеет состояние и соседей (это все ячейки, примыкающие к этой ячейке).
- Состояние каждой ячейки изменяется через определенное количество дискретных временных шагов в соответствии с набором правил.
- Состояние каждой ячейки на следующем временном шаге определяется текущим состоянием ячейки и состоянием соседних ячеек.
- Набор правил определяет обновление состояния ячейки.
Сейчас картина может выглядеть слишком абстрактной, поэтому попробуем представить проблему лесных пожаров с помощью КА. Для этого я буду использовать двумерную сетку. Каждая ячейка сетки будет принимать одно из следующих трех значений: дерево (1), лесной пожар (2), выжженная площадь (0). На рис. 1 показана случайно сгенерированная сетка 3×3, заполненная ячейками со всеми тремя состояниями.
Определим правила для нашей маленькой модели:
- Если возле ячейки с деревом есть лесной пожар хотя бы в одной из соседних ячеек, она начинает гореть.
- Если в ячейке есть дерево и рядом нет ни одного горящего соседа, то клетка начнет гореть с произвольной вероятностью f.
- Если ячейка пуста, в ней может вырасти дерево с произвольной вероятностью p.
Надеюсь, что рис. 1 и 2 убедили вас в том, что эта модель интуитивно понятна и проста, поэтому попробуем запрограммировать ее, используя только numpy и matplotlib.
Для этого проекта понадобится импортировать следующее:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib import colors
from matplotlib.patches import Rectangle
from datetime import timedelta,datetime
Поскольку мы работаем на Python, построим КА с помощью объектно-ориентированного программирования, то есть создадим класс. Пока что наш класс выглядит очень просто и выполняет всего три функции: инициализация ландшафта (сетки), проверка наличия пожара в соседних ячейках и распространение пожара:
class WFSim:
def __init__(self, f=0.01, p=1e-4, h=16, w=16):
self.f = f #вероятность пожара
self.p = p #вероятность вырастания дерева
self.h = h #высота сетки
self.w = w #ширина сетки
self.landscape = np.random.randint(1,2,(self.h,self.w)) #генерация сетки
def fire_neighbors_check(self, idx, jdx):
check = False
offsets = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
#offsets = координаты всех соседей ячейки (0,0)
for di,dj in offsets: #проверка, нет ли у кого-нибудь из соседей пожара
ni, nj = idx+di, jdx+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w:
if self.landscape[ni, nj] == 2:
check +=True
return check
def step(self, step):
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if self.landscape[i, j]==2:
new_landscape[i, j] = 0 #пепел после пожара
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #вырастание дерева
if (self.f > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #зажигание дерева
self.landscape = new_landscape.copy()
Теперь анимируем результаты с помощью matplotlib:
def update(frame):
im.set_data(Sim.landscape)
ax.axis('off')
Sim.step(frame)
return [im]
colors_list = ['black', 'forestgreen', 'orange'] #установка цветов ячейки
cmap = colors.ListedColormap(colors_list)
bounds = [0,1,2,3]
Sim = WFSim(h=64, w=64) #инициализация модели
fig, ax = plt.subplots(figsize=(16,16))
im = ax.imshow(Sim.landscape, cmap=cmap, norm=norm)
plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0,
hspace = 0, wspace = 0)
ani = animation.FuncAnimation(fig, update, frames=80, interval=20) #фреймы - количество шагов в модели
ani.save('simple.gif', fps=1.5, savefig_kwargs={'pad_inches':0})
plt.show()
Выглядит красиво, но проблема в том, что эта модель очень далека от реальности. Усложним ее!
Больше состояний ячейки
Первый шаг — добавить три новых состояния ячеек: вода, твердая порода, трава. Для удобства я изменю обозначения: вода (-3), твердая порода (-2), выжженная площадь (-1), трава (0), дерево (1), огонь (2). Все три новых класса будут инициализированы с самого начала, и в дальнейшем во время моделирования ячейки «вода» и «твердая порода» не смогут принимать другие значения (они фиксированы). Ячейки «трава» будут подчиняться следующим правилам:
- ячейка «трава» не может сгореть;
- сгоревшая ячейка может превратиться только в траву с вероятностью grass;
- ячейка «трава» может превратиться в ячейку «дерево» со (старой) вероятностью p (до этого мы превращали сгоревшую ячейку в дерево).
Таким образом, я попытался сделать модель более реалистичной. В природе после лесного пожара первым делом вырастает трава, только потом деревья, а вода и твердая порода не горят.
Итак, попробуем внедрить новые изменения в код:
class WFSim:
def __init__(self, f=0.01, p=1e-4, bedrock=0.005, water=0.05, h=16, w=16):
self.f = f #вероятность пожара
self.p = p #вероятность вырастания дерева
self.h = h #высота сетки
self.w = w #ширина сетки
self.bedrock = bedrock #вероятность твердой породы
self.water = water #вероятность воды
self.landscape = np.random.randint(0,2,(self.h,self.w)) #0 - трава, 1 - дерево
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
coef = 5 if self.surf_neighbors_check(i, j, "B") else 1
if self.bedrock*coef>np.random.random():
self.landscape[i, j] = -2
coef = 10 if self.surf_neighbors_check(i, j, "W") else 0.1
if self.water*coef>np.random.random():
self.landscape[i, j] = -3
def fire_neighbors_check(self, idx, jdx):
check = False
offsets = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
#offsets = координаты всех соседей ячейки (0,0)
for di,dj in offsets: #проверка, нет ли у кого-нибудь из соседей пожара
ni, nj = idx+di, jdx+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w:
if self.landscape[ni, nj] == 2:
check +=True
return check
def step(self, step):
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if self.landscape[i, j]==2:
new_landscape[i, j] = 0 #пепел после пожара
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #вырастание дерева
if (self.f > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #зажигание дерева
self.landscape = new_landscape.copy()
Как видите, процесс инициализации сетки изменился. На первом шаге мы случайным образом создаем сетку, состоящую только из травы и деревьев (соотношение «трава/деревья» нас не очень волнует). Затем, перебирая все ячейки с определенной вероятностью, мы присваиваем каждой ячейке значение воды или твердой породы. Но в природе скальные породы и водные бассейны часто существуют в кластерах. Поэтому, чтобы смоделировать большие скопления воды/пород, я ввел коэффициенты. Если у ячейки есть сосед с водой/твердой породой, то вероятность того, что в ней есть вода/твердая порода, увеличивается в 10 или 5 раз соответственно.
colors_list = ['steelblue', 'grey','black', 'olivedrab', 'forestgreen', 'orange']
cmap = colors.ListedColormap(colors_list)
bounds = [-3, -2, -1, 0,1,2,3]
norm = colors.BoundaryNorm(bounds, cmap.N)
Sim = WFSim(h=32, w=32)
plt.imshow(Sim.landscape, cmap=cmap, norm=norm)
plt.axis(False)
plt.savefig('Landscape.png', pad_inches=0, bbox_inches='tight')
plt.show()
Запустим то, что у нас получилось:
После добавления новых состояний ячеек можем пойти дальше и расширить внешние правила для этой симуляции. Конечно, я имею в виду погоду. Интенсивность лесных пожаров зависит от температуры воздуха, ветра и осадков. Так что попробуем реализовать эти три функции шаг за шагом и начнем с самой простой — температуры.
Температура
Температура сильно зависит от временной составляющей. Поэтому в дальнейшем я буду говорить о шагах моделирования как о часах. Первый шаг — это 00 часов по местному времени. Чтобы включить температуру в схему, нужно синтетически сгенерировать ее. Это несложно сделать, поскольку временной ряд температуры за один день имеет форму колокола (с максимальным значением около полудня). Затем можем использовать функцию синуса:
average_temp = 20
amplitude = 5
noise_level = 2
hours_in_day = 24
hours = np.arange(hours_in_day)
temperatures = average_temp + amplitude * np.sin(2 * np.pi * hours / 24 - np.pi / 2)
temperatures += np.random.normal(0, noise_level, 24) #шум
plt.plot(hours, temperatures, label='Simulated Temperature')
plt.xlabel('Hour of the Day')
plt.ylabel('Temperature (°C)')
plt.title('Simulated Average Daily Hourly Temperature')
plt.xticks(np.arange(0, 25, 3))
plt.legend()
plt.grid(True)
plt.savefig('temperature.png')
plt.show()
Таким образом, когда мы запускаем этот код, получаем 24 последовательных значения температуры воздуха. Как использовать их в нашей модели? Легко! Например, предположим, что если температура на определенном шаге превысит 25 °C, то вероятность f (начала горения какого-то дерева) удвоится. Заменим предыдущий код:
def step(self, step):
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if self.landscape[i, j]==2:
new_landscape[i, j] = -1 #пепел после пожара
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #вырастание дерева
if (self.f > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #зажигание дерева
self.landscape = new_landscape.copy()
На этот:
def step(self, step):
if step%24==0 and step>0:
self.temp = self.temperature()
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if self.landscape[i, j]==2:
new_landscape[i, j] = -1 #пепел после пожара
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #вырастание дерева
coef = 2 if self.temp[step%24]>25 else 1
if (self.f*coef > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #зажигание дерева
self.landscape = new_landscape.copy()
Ветер
Добавление ветра в эту модель потребует чуть больше усилий. В нашем случае ветер будет иметь только компонент направления (без скорости). Возможными направлениями ветра будут S (южный), N (северный), W (западный), E (восточный), NE (северо-восточный), NW (северо-западный), SE (юго-восточный), SW (юго-западный). В общем случае правило для ветра звучит так: огонь распространяется по ячейкам в соответствии с направлением ветра (рис. 8).
Чтобы реализовать это, нужно немного изменить правило смещения:
#....................
self.wind = np.random.choice(['calm', 'S', 'N', 'W', 'E', 'SW', 'SE', 'NW', 'NE'])
self.offsets = {'calm': [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)],
'N': [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1)],
'S': [(0,-1),(0,1),(1,-1),(1,0),(1,1)],
'E': [(-1,0),(-1,1),(0,1),(1,0),(1,1)],
'W': [(-1,-1),(-1,0),(0,-1),(1,-1),(1,0)],
'NE': [(-1,0),(-1,1),(0,1)],
'NW': [(-1,-1),(-1,0),(0,-1)],
'SE': [(0,1),(1,0),(1,1)],
'SW': [(0,-1),(1,-1),(1,0)]}
def fire_neighbors_check(self, idx, jdx):
check = False
offsets = self.offsets[self.wind]
for di,dj in offsets: #проверка, нет ли у соседей пожара
ni, nj = idx+di, jdx+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w:
if self.landscape[ni, nj] == 2:
check +=True
return check
Теперь можно выбрать направление ветра вручную или случайным образом и получить что-то вроде этого:
Осадки
Что является источником осадков? Правильно, облака. Поэтому введем в модель еще один тип ячеек: «облако (3)». Правило, связанное с облаками, будет звучать следующим образом:
Если облако накрывает сгоревшую или горящую ячейку, то на следующем шаге эта ячейка превращается в траву.
Теперь нам предстоит решить две сложные задачи:
- Реалистично сгенерировать облака.
- Перемещать облака.
Для генерации облаков введем новый параметр cloud, который представляет вероятность того, что на текущем шаге будет сгенерировано облако. Для генерации облака создадим новую функцию внутри класса. Идея будет заключаться в следующем:
- Случайным образом выбираем ячейку на ландшафте и присваиваем ей состояние cloud («облако» (3).
- Случайным образом выбираем одну соседнюю ячейку и присваиваем ей состояние «облако».
- Случайным образом выбираем одну соседнюю ячейку из ранее выбранных и также присваиваем ей состояние «облако».
- Повторите это N раз (N — размер облака).
Алгоритм позволяет создавать довольно реалистичные облака, которые имеют сложную форму (не просто квадраты или прямоугольники):
def generate_cloud(self):
size = 16 #произвольный размер облака
mask = np.zeros((self.h, self.w))
offsets = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
idx_prev, jdx_prev = np.random.randint(0,self.h),np.random.randint(0,self.w)
for _ in range(size):
for di,dj in [offsets[np.random.randint(len(offsets))]]:
ni, nj = idx_prev+di, jdx_prev+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w:
cell = (np.random.randint(min(0,idx_prev-1),min(self.h,idx_prev+1)),
np.random.randint(min(0,jdx_prev-1),min(self.w,jdx_prev+1)))
mask[ni, nj] = 1
idx_prev, jdx_prev = ni, nj
return mask.astype(bool)
Для перемещения облака потребуется ввести новую переменную, сохраняющую ландшафт на предыдущем этапе. Зачем это нужно? Представьте облако, накрывающее ячейку с водой, ячейку с деревом и сожженную ячейку. Когда оно перемещается, нужно оставлять ячейки с водой и деревьями прежними, но изменять состояние сожженной ячейки на состояние «трава». Поэтому сохранение предыдущего состояния неизбежно.
Облака также будут двигаться по направлению ветра, поэтому на каждой итерации облако перемещается с шагом в одну ячейку в соответствии с направлением ветра.
Если ветра нет, облака распадаются и через один шаг исчезают.
Напишем для этого код!
#..........
self.cloud_offsets = {'calm': [],
'N': [(1,0)],
'S': [(-1,0)],
'E': [(0,-1)],
'W': [(0,1)],
'NE': [(1,-1)],
'NW': [(1,1)],
'SE': [(-1,-1)],
'SW': [(-1,1)]}
def cloud_move(self):
offsets = self.cloud_offsets[self.wind]
mask = np.zeros((self.h, self.w))
for idx in range(self.landscape.shape[0]):
for jdx in range(self.landscape.shape[1]):
for di,dj in offsets:
ni, nj = idx+di, jdx+dj
if nj>=0 and ni>=0 and ni<self.h and nj<self.w and self.landscape[idx, jdx]==3:
mask[ni,nj] = 1
return mask.astype(bool)
Теперь можем добавить в функцию step генерацию и движение облаков:
def step(self, step):
if step%24==0 and step>0:
self.temp = self.temperature()
new_landscape = self.landscape.copy()
for i in range(self.landscape.shape[0]):
for j in range(self.landscape.shape[1]):
if new_landscape[i, j] == 3:
if self.old_landscape[i,j]==-1 or self.old_landscape[i,j]==2:
new_landscape[i, j] = 0 #вырастание травы, если ячейка сгорела или горит
else:
new_landscape[i, j] = self.old_landscape[i,j] #иначе пропустить
if new_landscape[i, j] == 2:
new_landscape[i, j] = -1 #пепел
if self.p > np.random.rand() and self.landscape[i, j]==0:
new_landscape[i, j] = 1 #вырастить дерево в ячейке "трава"
coef = 2 if self.temp[step%24]>25 else 1
if (self.f*coef > np.random.rand() or self.fire_neighbors_check(i, j)) and self.landscape[i, j]==1:
new_landscape[i, j] = 2 #горение
self.old_landscape = new_landscape.copy() #сохранить обновленный ландшафт
if 3 in self.landscape and self.wind!='calm':
new_landscape[self.cloud_move()] = 3 #перемещать облако, если оно есть
if (self.cloud > np.random.rand()):
new_landscape[self.generate_cloud()] = 3 #создать новое облако
self.landscape = new_landscape.copy() #обновить основной ландшафт
Метрики
Мне очень нравится то, что у нас получилось, но осталось добавить последнюю деталь — цифры. Создадим небольшое поле, показывающее температуру, ветер, наличие деревьев, соотношение сгоревших пикселей и шаг:
def update(frame):
im.set_data(Sim.landscape)
info_text = (
f"Date: {initial_date + timedelta(hours=frame)}\n"
f"Wind: {Sim.wind}\n"
f"Temperature: {round(Sim.temp[frame%24],1)}°C\n"
f"Burned area: {Sim.burned_ratio} %\n"
f"Tree cover: {Sim.tree_cover} %"
)
text.set_text(info_text)
ax.axis('off')
Sim.step(frame)
return [im]
Sim = WFSim(h=64, w=64)
initial_date = datetime(2012, 8, 19, 0)
fig, ax = plt.subplots(figsize=(16,16))
im = ax.imshow(Sim.landscape, cmap=cmap, norm=norm)
text = ax.text(2, 8, "", ha='left', fontsize=20, color='white',
bbox=dict(facecolor='black', alpha=0.8, edgecolor='black',boxstyle='round,pad=0.5'))
plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0,
hspace = 0, wspace = 0)
ani = animation.FuncAnimation(fig, update, frames=80, interval=20)
ani.save('final.gif', fps=1.5, savefig_kwargs={'pad_inches':0})
plt.show()
Теперь намного лучше!
Эта модель уже сложнее, но ее, конечно, можно усовершенствовать, введя новые внешние факторы/состояния ячеек, так что она будет еще более приближенной к условиям реального мира.
Читайте также:
- Визуализация параметров градиентного спуска в Torch
- Лучшая IDE для Python-разработки в 2024 году
- Битва 4 инструментов визуализации данных на языке Python
Читайте нас в Telegram, VK и Дзен
Перевод статьи Aleksei Rozanov: Simulating Wildfires in a Forest
