使用 Python 实现随机中点位移法生成逼真的裂隙面

一、随机中点位移法简介
- 1. 什么是随机中点位移法？
- 2. 应用领域
二、 Python 代码实现
- 1. 导入必要的库
- 2. 函数定义：随机中点位移法核心逻辑
- 3. 设置随机数种子
- 4. 初始化二维裂隙面
- 5. 初始化网格的四个顶点
- 6. 初始化步长
- 7. 开始迭代生成裂隙面
- 8. 钻石步（ Diamond Step ）
- 9. 方形步（ Square Step ）
- 10. 更新步长和粗糙度
- 11. 返回最终结果
- 12. 调用函数并可视化结果
三、完整代码
四、结果展示与分析
五、扩展与应用
- 1. 三维裂隙面生成
- 2. 工程领域应用
- 3. 性能优化

裂隙面在地质学、岩石力学、计算机图形学等领域有着广泛的应用。它们的随机性和分形特性使得模拟真实裂隙面成为一项重要的研究课题。本文将通过 Python ，结合经典的随机中点位移法（ Random Midpoint Displacement , RMD ），逐步创建一个二维的随机裂隙面，并详细解析每一步代码的实现过程。

先看下效果图：
在这里插入图片描述

一、随机中点位移法简介

1. 什么是随机中点位移法？

随机中点位移法是一种简单且高效的分形生成算法。它的基本思想是从一个大尺度开始，通过不断细化网格，并在每次细化时为中间点添加随机扰动，从而生成具有大尺度结构和小尺度细节的表面。

这一方法的特点是：

高效性：算法复杂度较低，适合生成大规模的表面。
随机性：通过随机扰动，模拟自然界中裂隙、地形等不规则表面的特性。
分形特性：生成的结果具有分形维数，能很好地模拟自然界的粗糙表面。

2. 应用领域

地质学：生成裂隙面，用于研究地下水流动、岩石的强度和稳定性。
计算机图形学：生成自然地貌（如山脉、地形）或粗糙材质。
工程应用：模拟裂隙网络，分析裂隙对渗透率的影响。

二、 Python 代码实现

让我们直接进入代码部分，并通过逐行详解，帮助你全面掌握随机中点位移法的实现。

1. 导入必要的库

import numpy as np
import matplotlib.pyplot as plt

numpy ：用于处理多维数组和数学计算。我们将借助它来创建二维数组并进行随机数生成。
matplotlib.pyplot ：用于可视化生成的裂隙面。这是 Python 中非常强大的绘图库。

2. 函数定义：随机中点位移法核心逻辑

def midpoint_displacement_2d(size, roughness, seed=None):"""使用随机中点位移法生成二维粗糙裂隙面。参数:- size: 裂隙面的网格尺寸（必须是2的幂次方+1，例如65, 129, 257等）。- roughness: 粗糙度参数，控制表面的起伏程度。- seed: 随机数种子，用于生成可重复的结果。返回值:- surface: 生成的二维裂隙面。"""

函数说明

size ：二维数组的大小，必须为 2^n + 1 （如 129 , 257 等）。这种结构使网格能够在每次迭代中被整齐地细分。
roughness ：粗糙度参数，控制生成裂隙面的起伏程度。值越大，表面越粗糙；值越小，表面越平滑。
seed ：随机数种子，用于确保生成的结果是可重复的。指定种子后，每次运行都会生成相同的裂隙面。

3. 设置随机数种子

    if seed is not None:np.random.seed(seed)

解释

如果用户提供了 seed 值，那么我们使用 np.random.seed(seed) 设置随机数生成器的种子。
这样可以确保生成的随机数序列是固定的，从而保证结果的可重复性。

4. 初始化二维裂隙面

    surface = np.zeros((size, size))

np.zeros((size, size)) ：创建一个大小为 (size, size) 的二维数组，初始值全为 0 。
这个数组将表示裂隙面的高度值，每一个元素对应网格的一个节点。

5. 初始化网格的四个顶点

    surface[0, 0] = np.random.uniform(-1, 1)surface[0, -1] = np.random.uniform(-1, 1)surface[-1, 0] = np.random.uniform(-1, 1)surface[-1, -1] = np.random.uniform(-1, 1)

解释

随机为二维数组的四个角点赋值，值范围为 [-1, 1] 。
这四个顶点值将作为整个裂隙面生成的初始条件。

6. 初始化步长

    step_size = size - 1

解释

设置初始步长为 size - 1 。
步长是指当前迭代中，用于分割网格的大小。在每次迭代中，步长会减半，从而逐步细化网格。

7. 开始迭代生成裂隙面

    while step_size > 1:half_step = step_size // 2

解释

使用 while 循环，直到步长小于或等于 1 。
在每次循环中，计算半步长 half_step = step_size // 2 ，用于在网格的中点和边界点之间插值。

8. 钻石步（ Diamond Step ）

        for x in range(0, size - 1, step_size):for y in range(0, size - 1, step_size):avg = (surface[x, y]+ surface[x + step_size, y]+ surface[x, y + step_size]+ surface[x + step_size, y + step_size]) / 4surface[x + half_step, y + half_step] = avg + np.random.uniform(-1, 1) * roughness

说明

目标：为每个网格单元的中心点生成一个新值。
avg ：计算当前网格四个顶点的平均值。
np.random.uniform(-1, 1) * roughness ：添加一个随机扰动，其幅度由 roughness 控制。
更新后的值存储在网格的中心点位置 surface[x + half_step, y + half_step] 。

9. 方形步（ Square Step ）

        for x in range(0, size - 1, half_step):for y in range((x + half_step) % step_size, size - 1, step_size):avg = 0count = 0if x - half_step >= 0:avg += surface[x - half_step, y]count += 1if x + half_step < size:avg += surface[x + half_step, y]count += 1if y - half_step >= 0:avg += surface[x, y - half_step]count += 1if y + half_step < size:avg += surface[x, y + half_step]count += 1surface[x, y] = avg / count + np.random.uniform(-1, 1) * roughness

说明

目标：为网格的每条边界点生成一个新值。
通过计算边界点相邻的点的平均值，并添加随机扰动，为边界点赋值。

10. 更新步长和粗糙度

        step_size //= 2roughness /= 2

说明

每次迭代后，步长减半，逐步细化网格。
同时，降低随机扰动幅度（粗糙度），使得细节更平滑。

11. 返回最终结果

    return surface

说明

返回生成的裂隙面（二维数组）。

12. 调用函数并可视化结果

size = 129
roughness = 1.0
seed = 42
surface = midpoint_displacement_2d(size, roughness, seed)plt.figure(figsize=(10, 8))
plt.imshow(surface, cmap='terrain', origin='upper')
plt.colorbar(label="Height")
plt.title("Random Midpoint Displacement Surface")
plt.xlabel("X-axis")
plt.ylabel("Y-axis")
plt.show()

说明

参数：
- size=129 ：网格大小为 129 × 129 。
- roughness=1.0 ：初始粗糙度。
- seed=42 ：随机种子，用于生成可复现的结果。
使用 matplotlib 绘制结果：
- imshow ：展示二维裂隙面。
- cmap='terrain' ：用地形色彩映射模拟裂隙表面。

三、完整代码

import numpy as np
import matplotlib.pyplot as pltdef midpoint_displacement_2d(size, roughness, seed=None):"""使用随机中点位移法生成二维粗糙裂隙面。参数:- size: 裂隙面的网格尺寸（必须是2的幂次方+1，例如65, 129, 257等）。- roughness: 粗糙度参数，控制表面的起伏程度。- seed: 随机数种子，用于生成可重复的结果。返回值:- surface: 生成的二维裂隙面。"""if seed is not None:np.random.seed(seed)# 初始化一个二维数组，并将四个顶点随机赋值surface = np.zeros((size, size))surface[0, 0] = np.random.uniform(-1, 1)surface[0, -1] = np.random.uniform(-1, 1)surface[-1, 0] = np.random.uniform(-1, 1)surface[-1, -1] = np.random.uniform(-1, 1)step_size = size - 1  # 初始步长while step_size > 1:half_step = step_size // 2# 钻石步：计算中点值for x in range(0, size - 1, step_size):for y in range(0, size - 1, step_size):avg = (surface[x, y]+ surface[x + step_size, y]+ surface[x, y + step_size]+ surface[x + step_size, y + step_size]) / 4surface[x + half_step, y + half_step] = avg + np.random.uniform(-1, 1) * roughness# 方形步：计算边点值for x in range(0, size - 1, half_step):for y in range((x + half_step) % step_size, size - 1, step_size):avg = 0count = 0if x - half_step >= 0:avg += surface[x - half_step, y]count += 1if x + half_step < size:avg += surface[x + half_step, y]count += 1if y - half_step >= 0:avg += surface[x, y - half_step]count += 1if y + half_step < size:avg += surface[x, y + half_step]count += 1surface[x, y] = avg / count + np.random.uniform(-1, 1) * roughness# 缩小步长，并降低粗糙度比例step_size //= 2roughness /= 2return surface# 参数设置
size = 129  # 网格尺寸（2^n + 1，例如65, 129, 257）
roughness = 1.0  # 粗糙度系数
seed = 42  # 随机种子（可选）# 生成裂隙面
surface = midpoint_displacement_2d(size, roughness, seed)# 绘制裂隙面
plt.figure(figsize=(10, 8))
plt.imshow(surface, cmap='terrain', origin='upper')
plt.colorbar(label="Height")
plt.title("Random Midpoint Displacement Surface")
plt.xlabel("X-axis")
plt.ylabel("Y-axis")
plt.show()