识别3D NumPy孔隙度数组中贯穿各层同类孔隙区的首层网格单元

一、先解决“统计层数为0”的核心问题 #

你遇到的统计结果全为0，基本是两个原因：要么分类阈值和孔隙度数据范围不匹配，要么分类/统计的轴方向搞反了。下面是修正后的完整代码，先把基础逻辑跑通：

import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt

# 1. 生成符合储层特征的孔隙度数据（0-1范围，模拟真实分布）
porosity_array = np.random.normal(loc=0.3, scale=0.15, size=(40, 50, 50))
porosity_array = np.clip(porosity_array, 0.1, 0.8)  # 限制合理范围，避免极端值

# 2. 高斯平滑处理
smoothed_poro = gaussian_filter(porosity_array, sigma=(1, 2, 2))  # 层间平滑弱，平面内平滑强

# 3. 用分位数设置分类阈值（比固定值更适配数据分布）
low_thresh = np.percentile(smoothed_poro, 33)
mid_thresh = np.percentile(smoothed_poro, 67)

# 4. 分类：(层,行,列) → 每个网格点的类别（0=差，1=中，2=优）
poro_classes = np.zeros_like(smoothed_poro, dtype=int)
poro_classes[smoothed_poro >= mid_thresh] = 2  # 优（高孔隙度）
poro_classes[(smoothed_poro >= low_thresh) & (smoothed_poro < mid_thresh)] = 1  # 中
poro_classes[smoothed_poro < low_thresh] = 0  # 差（低孔隙度）

# 5. 统计每个首层网格（第0层）垂直方向穿过各类别的层数
# 注意：axis=0是层方向，统计每个(行,列)点在40层中属于各类别的次数
count_low = np.sum(poro_classes[:, :, :] == 0, axis=0)
count_mid = np.sum(poro_classes[:, :, :] == 1, axis=0)
count_high = np.sum(poro_classes[:, :, :] == 2, axis=0)

# 验证：打印首层几个点的统计结果，应该不会全为0
print("首层随机点的低孔隙度层数：", count_low[10,10])
print("首层随机点的中孔隙度层数：", count_mid[10,10])
print("首层随机点的高孔隙度层数：", count_high[10,10])

二、高效实现井位识别（分散+垂直穿透多） #

要同时满足“垂直穿透对应类别多”和“井位不聚集”，可以用贪心筛选+空间约束的思路：

步骤1：对每个类别生成优先级图 #

比如针对高孔隙度（优）类别，直接拿count_high作为优先级，值越高的点越优先选。

步骤2：贪心选点+排除邻域 #

def select_wells(priority_map, num_wells, exclude_radius=3):
    """
    从优先级图中选分散的井位
    :param priority_map: (50,50)的优先级数组（值越高越优先）
    :param num_wells: 要选的井位数
    :param exclude_radius: 选点后排除周围N个网格
    :return: 选中的井位坐标列表[(行,列), ...]
    """
    wells = []
    priority_copy = priority_map.copy()
    
    for _ in range(num_wells):
        # 找当前优先级最高的点
        max_idx = np.unravel_index(np.argmax(priority_copy), priority_copy.shape)
        wells.append(max_idx)
        
        # 排除该点周围radius范围内的所有点，避免聚集
        r = exclude_radius
        row, col = max_idx
        # 限制边界，防止越界
        row_min = max(0, row - r)
        row_max = min(priority_copy.shape[0], row + r + 1)
        col_min = max(0, col - r)
        col_max = min(priority_copy.shape[1], col + r + 1)
        priority_copy[row_min:row_max, col_min:col_max] = -np.inf  # 设为极小值，不再选中
    
    return wells

# 示例：选5个高孔隙度优先的分散井位
high_wells = select_wells(count_high, num_wells=5, exclude_radius=4)
print("选中的高孔隙度井位：", high_wells)

三、Jet色图可视化优化 #

直接用matplotlib的jet色图对应分类，同时叠加井位标记：

# 1. 首层孔隙度分类可视化
plt.figure(figsize=(8,8))
im = plt.imshow(poro_classes[0, :, :], cmap='jet', vmin=0, vmax=2)
# 标记选中的井位
for (row, col) in high_wells:
    plt.scatter(col, row, marker='*', s=200, c='white', edgecolor='black')
plt.colorbar(im, ticks=[0,1,2], label='孔隙度类别：0=差(蓝),1=中(黄),2=优(红)')
plt.title('首层孔隙度分类与井位分布')
plt.show()

# 2. 垂直剖面可视化（比如沿着某口井的列方向）
well_row, well_col = high_wells[0]
plt.figure(figsize=(10,4))
im = plt.imshow(poro_classes[:, well_row, :], cmap='jet', aspect='auto', vmin=0, vmax=2)
plt.scatter(well_col, 0, marker='*', s=200, c='white', edgecolor='black')
plt.colorbar(im, ticks=[0,1,2], label='孔隙度类别')
plt.title(f'井位({well_row},{well_col})垂直剖面孔隙度分布')
plt.xlabel('列')
plt.ylabel('层')
plt.show()

关键优化点总结 #

• 用分位数阈值替代固定值：避免因随机数据分布变化导致分类失效，更适配真实储层数据的统计特征
• 向量化统计：完全用numpy的sum(axis=0)替代循环，效率提升几个数量级
• 贪心选点+邻域排除：简单高效，保证井位分散，比聚类算法更容易理解和调整
• 可视化叠加标记：直接验证井位是否落在目标区域，垂直剖面可直观看到穿透层数

内容的提问来源于stack exchange，提问作者Skylimit