CycleUser1. 核心概念：临界与混沌
通常我们认为，一个复杂的系统（如森林大火）必然需要复杂的机制来解释。但物理学家发现，哪怕是最简单的规则，也能涌现出极其复杂的行为。这在森林火灾模型（Forest Fire Model）中体现得淋漓尽致。
回归到最自然的森林生态循环：树木生长、雷电引燃、火势蔓延，这个模型展示了一个非常重要的概念：自组织临界性（Self-Organized Criticality, SOC）。系统会自动演化到一个"临界状态"，在这个状态下，哪怕是一个小小的火星（微小扰动），既可能瞬间熄灭，也可能引发一场烧毁半个森林的超级大火（雪崩效应）。
1.1 角色与规则
我们把森林看作一个二维网格，每个格子有三种状态：
- 空地（Empty, 0）：这里什么也没有（可能是刚烧完留下的灰烬）。
- 树木（Tree, 1）：一棵健康的树。
- 燃烧（Burning, 2）：正在燃烧的树。
1.2 动力学演化
时间一步步向前推进，每一步，每个格子都会根据以下规则同时更新：
- 树木生长：如果一个格子是空地，它有概率 $p$ 长出一棵新树。
- 火势蔓延：如果一棵树的上下左右四个邻居中，至少有一个在燃烧，这棵树就会被点燃。
- 雷电引燃：即使周围没有火，一棵树也有极小的概率 $f$ 被雷电击中而自燃。
- 燃尽熄灭：正在燃烧的树，在下一步会变成空地。
2. 数学模型：元胞自动机
这属于元胞自动机（Cellular Automata, CA）模型。与微分方程不同，它是离散的。
设 $S_{i,j}^{(t)}$ 为时间 $t$ 时位置 $(i,j)$ 的状态。
-
空地 $\to$ 树木: $$ S_{i,j}^{(t+1)} = 1 \quad \text{if } S_{i,j}^{(t)} = 0 \text{ with probability } p $$
-
树木 $\to$ 燃烧: $$ S_{i,j}^{(t+1)} = 2 \quad \text{if } S_{i,j}^{(t)} = 1 \text{ AND } (\exists \text{ neighbor } = 2 \text{ OR with probability } f) $$
-
燃烧 $\to$ 空地: $$ S_{i,j}^{(t+1)} = 0 \quad \text{if } S_{i,j}^{(t)} = 2 $$
关键参数比率
有趣的事情发生在树木生长速率 $p$ 远大于雷电频率 $f$ 时（$f \ll p$）。 * 树木慢慢生长，森林越来越茂密，形成连通的"燃料簇"。 * 一旦雷电击中某棵树，火势会顺着连通的树木迅速蔓延。 * 如果森林不够密，火烧一小片就停了。 * 如果森林太密，一场大火会烧掉绝大部分，系统"重置"。 * 系统会在这两者之间摆动，长期处于一种临界状态。
3. Python 代码演示
我们将使用 numpy 进行高效的网格计算，并使用 PySide6 (Qt for Python) 构建一个功能完备的交互式仿真程序。该程序完全复刻了经典的 Drossel-Schwabl 模型特性，并增加了高级控制功能：
- 参数实时调节：精确控制树木生长率 $p$ (0.00-0.10) 和雷电频率 $f$ (1e-6量级)。
- 自适应变速 (Adaptive Speed)：这是模拟 SOC 现象的关键。程序可以设置"平时"（只有生长）快进运行，一旦检测到火灾（Fire），自动降速到逐帧播放。这让你既不需要等待漫长的树木生长，又不会错过精彩的火灾蔓延细节。
- 动态录制：支持将模拟过程录制并导出为 GIF 动画。
图示：黑色为空地，绿色为树木，橙色为正在燃烧的火线。
完整代码
你需要安装 numpy, PySide6, Pillow 库：
pip install numpy PySide6 Pillow
import sys
import numpy as np
from PIL import Image
from PySide6.QtWidgets import (QApplication, QMainWindow, QWidget, QVBoxLayout,
QHBoxLayout, QPushButton, QLabel, QSlider,
QSpinBox, QFileDialog, QMessageBox, QGroupBox, QCheckBox)
from PySide6.QtCore import QTimer, Qt
from PySide6.QtGui import QImage, QPainter
# --- Model Logic ---
class ForestFireModel:
def __init__(self, width=200, height=200, p=0.01, f=0.00001, seed=None):
self.width = width
self.height = height
self.p = p
self.f = f
self.seed = seed
# Use a dedicated random number generator
self.rng = np.random.default_rng(seed)
self.grid = np.zeros((height, width), dtype=int)
self.reset()
def reset(self):
if self.seed is not None:
self.rng = np.random.default_rng(self.seed)
# Initialize with some random trees (50% density)
self.grid = self.rng.choice([0, 1], size=(self.height, self.width), p=[0.5, 0.5])
def step(self):
# Create a copy to update synchronously
new_grid = self.grid.copy()
# Masks for current states
empty_mask = (self.grid == 0)
tree_mask = (self.grid == 1)
fire_mask = (self.grid == 2)
# Rule 1: Growth (Empty -> Tree)
growth_random = self.rng.random((self.height, self.width))
new_grid[empty_mask & (growth_random < self.p)] = 1
# Rule 2: Fire Spread (Tree -> Fire if neighbor burning)
# Check 4-neighbor (von Neumann)
north = np.roll(fire_mask, 1, axis=0)
south = np.roll(fire_mask, -1, axis=0)
west = np.roll(fire_mask, 1, axis=1)
east = np.roll(fire_mask, -1, axis=1)
neighbor_on_fire = (north | south | west | east)
# Rule 3: Lightning (Tree -> Fire random)
lightning_random = self.rng.random((self.height, self.width))
lightning_strike = (lightning_random < self.f)
# Trees ignite if neighbor burns OR lightning strikes
ignite_mask = tree_mask & (neighbor_on_fire | lightning_strike)
new_grid[ignite_mask] = 2
# Rule 4: Burnt (Fire -> Empty)
new_grid[fire_mask] = 0
self.grid = new_grid
# Return True if any fire exists
return np.any(new_grid == 2)
# --- Custom Widget for Visualization ---
class SimulationCanvas(QWidget):
def __init__(self, model):
super().__init__()
self.model = model
# Color Map: 0=Black, 1=Green, 2=OrangeRed
self.colors = np.array([
[0, 0, 0], # Empty
[34, 139, 34], # Tree
[255, 69, 0] # Fire
], dtype=np.uint8)
self.image = None
def update_canvas(self):
# Map grid values to RGB colors
rgb_grid = self.colors[self.model.grid]
# Create QImage from the RGB data
height, width, channel = rgb_grid.shape
bytes_per_line = 3 * width
# QImage references the data, so we must make a deep copy
temp_image = QImage(rgb_grid.data, width, height, bytes_per_line, QImage.Format_RGB888)
self.image = temp_image.copy()
self.update() # Trigger paintEvent
def paintEvent(self, event):
if self.image:
painter = QPainter(self)
# Scale image to fit the widget, keeping aspect ratio
rect = self.rect()
scaled_image = self.image.scaled(rect.size(), Qt.KeepAspectRatio, Qt.FastTransformation)
# Center the image
x = (rect.width() - scaled_image.width()) // 2
y = (rect.height() - scaled_image.height()) // 2
painter.drawImage(x, y, scaled_image)
def get_current_frame(self):
"""Returns the current grid as a PIL Image for export"""
if self.image:
rgb_grid = self.colors[self.model.grid]
return Image.fromarray(rgb_grid)
return None
# --- Main Window ---
class MainWindow(QMainWindow):
def __init__(self):
super().__init__()
self.setWindowTitle("Forest Fire Simulation (Drossel-Schwabl)")
self.resize(1100, 750)
# Default Parameters (Veritasium style)
self.default_p = 0.01
self.default_f = 0.00001
self.default_size = 200
self.default_seed = 42
self.recording = False
self.frames = []
self.steps_per_draw = 1
self.adaptive_speed = True
# Model
self.model = ForestFireModel(self.default_size, self.default_size, self.default_p, self.default_f, self.default_seed)
# Timer
self.timer = QTimer()
self.timer.timeout.connect(self.update_simulation)
self.timer.setInterval(10) # Default fast refresh
# UI Layout
main_widget = QWidget()
self.setCentralWidget(main_widget)
layout = QHBoxLayout(main_widget)
# Left: Controls
control_panel = QWidget()
control_layout = QVBoxLayout(control_panel)
control_panel.setFixedWidth(320)
layout.addWidget(control_panel)
# Right: Canvas
self.canvas = SimulationCanvas(self.model)
layout.addWidget(self.canvas, stretch=1)
# --- Controls ---
grp_sim = QGroupBox("Control")
l_sim = QVBoxLayout()
self.btn_start = QPushButton("Start")
self.btn_start.clicked.connect(self.toggle_simulation)
l_sim.addWidget(self.btn_start)
self.btn_reset = QPushButton("Reset")
self.btn_reset.clicked.connect(self.reset_simulation)
l_sim.addWidget(self.btn_reset)
# Seed Control
l_seed = QHBoxLayout()
l_seed.addWidget(QLabel("Seed:"))
self.spin_seed = QSpinBox()
self.spin_seed.setRange(0, 999999)
self.spin_seed.setValue(self.default_seed)
self.spin_seed.valueChanged.connect(self.update_seed)
l_seed.addWidget(self.spin_seed)
l_sim.addLayout(l_seed)
grp_sim.setLayout(l_sim)
control_layout.addWidget(grp_sim)
# Group: Parameters
grp_params = QGroupBox("Parameters")
l_params = QVBoxLayout()
# Growth Rate (p)
l_params.addWidget(QLabel("Tree Growth Rate (p):"))
self.slider_p = QSlider(Qt.Horizontal)
self.slider_p.setRange(0, 1000) # 0.0 to 0.1
self.slider_p.setValue(int(self.default_p * 10000))
self.slider_p.valueChanged.connect(self.update_params)
l_params.addWidget(self.slider_p)
self.lbl_p = QLabel(f"{self.default_p:.4f}")
l_params.addWidget(self.lbl_p)
# Lightning Rate (f)
l_params.addWidget(QLabel("Lightning Rate (f):"))
self.slider_f = QSlider(Qt.Horizontal)
self.slider_f.setRange(0, 1000) # 0.0 to 0.001
self.slider_f.setValue(int(self.default_f * 1000000)) # Scale for 1e-6
self.slider_f.valueChanged.connect(self.update_params)
l_params.addWidget(self.slider_f)
self.lbl_f = QLabel(f"{self.default_f:.6f}")
l_params.addWidget(self.lbl_f)
# Grid Size
l_params.addWidget(QLabel("Grid Size (Square):"))
self.spin_size = QSpinBox()
self.spin_size.setRange(50, 500)
self.spin_size.setValue(self.default_size)
self.spin_size.setSingleStep(50)
self.spin_size.valueChanged.connect(self.resize_grid)
l_params.addWidget(self.spin_size)
grp_params.setLayout(l_params)
control_layout.addWidget(grp_params)
# Group: Speed Control
grp_speed = QGroupBox("Speed Control")
l_speed = QVBoxLayout()
# Refresh Interval
l_speed.addWidget(QLabel("Refresh Interval (ms):"))
self.slider_interval = QSlider(Qt.Horizontal)
self.slider_interval.setRange(1, 200)
self.slider_interval.setValue(10)
self.slider_interval.valueChanged.connect(self.update_speed_settings)
l_speed.addWidget(self.slider_interval)
self.lbl_interval = QLabel("10 ms")
l_speed.addWidget(self.lbl_interval)
# Steps per Draw
l_speed.addWidget(QLabel("Steps per Frame (Turbo):"))
self.slider_steps = QSlider(Qt.Horizontal)
self.slider_steps.setRange(1, 1000) # Allow up to 1000 steps per frame for high p
self.slider_steps.setValue(1)
self.slider_steps.valueChanged.connect(self.update_speed_settings)
l_speed.addWidget(self.slider_steps)
self.lbl_steps = QLabel("1 step/frame")
l_speed.addWidget(self.lbl_steps)
# Adaptive Speed
self.chk_adaptive = QCheckBox("Adaptive Speed (Slow for Fire)")
self.chk_adaptive.setChecked(True)
self.chk_adaptive.toggled.connect(self.update_speed_settings)
l_speed.addWidget(self.chk_adaptive)
grp_speed.setLayout(l_speed)
control_layout.addWidget(grp_speed)
# Group: Export
grp_export = QGroupBox("Export")
l_export = QVBoxLayout()
self.btn_record = QPushButton("Start Recording")
self.btn_record.setCheckable(True)
self.btn_record.clicked.connect(self.toggle_recording)
l_export.addWidget(self.btn_record)
self.lbl_record_status = QLabel("Frames: 0")
l_export.addWidget(self.lbl_record_status)
grp_export.setLayout(l_export)
control_layout.addWidget(grp_export)
control_layout.addStretch()
self.canvas.update_canvas()
def toggle_simulation(self):
if self.timer.isActive():
self.timer.stop()
self.btn_start.setText("Start")
else:
self.timer.start()
self.btn_start.setText("Pause")
def reset_simulation(self):
self.model.reset()
self.canvas.update_canvas()
def update_seed(self):
self.model.seed = self.spin_seed.value()
self.reset_simulation()
def update_params(self):
p_val = self.slider_p.value() / 10000.0
self.model.p = p_val
self.lbl_p.setText(f"{p_val:.4f}")
f_val = self.slider_f.value() / 1000000.0
self.model.f = f_val
self.lbl_f.setText(f"{f_val:.6f}")
def update_speed_settings(self):
interval = self.slider_interval.value()
self.timer.setInterval(interval)
self.lbl_interval.setText(f"{interval} ms")
self.steps_per_draw = self.slider_steps.value()
self.lbl_steps.setText(f"{self.steps_per_draw} step/frame")
self.adaptive_speed = self.chk_adaptive.isChecked()
def resize_grid(self):
size = self.spin_size.value()
self.model = ForestFireModel(size, size, self.model.p, self.model.f, self.model.seed)
self.canvas.model = self.model
self.canvas.update_canvas()
def toggle_recording(self):
if self.btn_record.isChecked():
self.recording = True
self.frames = []
self.btn_record.setText("Stop & Save GIF")
self.lbl_record_status.setText("Recording...")
else:
self.recording = False
self.btn_record.setText("Start Recording")
self.save_gif()
def save_gif(self):
if not self.frames: return
filename, _ = QFileDialog.getSaveFileName(self, "Save GIF", "", "GIF Files (*.gif)")
if filename:
if not filename.endswith('.gif'): filename += '.gif'
self.frames[0].save(filename, save_all=True, append_images=self.frames[1:], duration=50, loop=0)
QMessageBox.information(self, "Success", f"Saved to {filename}")
self.frames = []
self.lbl_record_status.setText("Frames: 0")
def update_simulation(self):
# Determine how many steps to run
steps_to_run = self.steps_per_draw
# Check for fire in current state to decide speed
has_fire = np.any(self.model.grid == 2)
if self.adaptive_speed and has_fire:
steps_to_run = 1 # Slow down to watch the fire
# Run steps
fire_still_exists = has_fire
for _ in range(steps_to_run):
fire_still_exists = self.model.step()
if self.adaptive_speed and fire_still_exists and steps_to_run > 1:
break
self.canvas.update_canvas()
if self.recording:
frame = self.canvas.get_current_frame()
if frame:
self.frames.append(frame)
self.lbl_record_status.setText(f"Frames: {len(self.frames)}")
if __name__ == "__main__":
app = QApplication(sys.argv)
window = MainWindow()
window.show()
sys.exit(app.exec())
4. 深入思考：为什么这很重要？
这个简单的模型揭示了自然界中许多灾难性事件的本质。
- 无法预测的大灾难：在临界状态下，火灾的大小服从幂律分布（Power Law）。这意味着，虽然大多数火灾都很小，但超级大火是不可避免的，而且你无法预测下一次雷击是会只烧掉一棵树，还是会烧掉整片森林。
- 生态平衡：火灾本身也是生态系统的一部分。它清理了老化的森林，为新树的生长腾出了空间（空地）。如果没有火灾，森林会变得过度拥挤，一旦发生火灾就是毁灭性的。
- 应用广泛：同样的原理可以用来解释地震（地壳应力积累与释放）、股市崩盘（市场情绪积累与恐慌抛售）甚至大脑神经元放电。
5. 参考文献 (References)
- Drossel, B., & Schwabl, F. (1992). Self-organized critical forest-fire model. Physical Review Letters, 69(11), 1629.
- 提出了这个经典的自组织临界森林火灾模型。
- Bak, P., Tang, C., & Wiesenfeld, K. (1987). Self-organized criticality: An explanation of the 1/f noise. Physical Review Letters, 59(4), 381.
- 关于自组织临界性（沙堆模型）的奠基性论文。
- Turcotte, D. L. (1999). Self-organized criticality. Reports on progress in physics, 62(10), 1377.
- 综述了 SOC 在地球物理学和自然灾害中的应用。