3D储层网格连通性与迂曲度计算代码运行过慢求助

核心性能瓶颈与逻辑问题 #

1. 迂曲度计算的性能与逻辑问题 #

• 使用NetworkX构建79万+节点的大型3D图，Python级别的三重循环添加边，效率极低。
• nx.all_pairs_shortest_path_length针对超大规模图的时间复杂度爆炸，完全无法完成计算。
• 重复计算渗透率与孔隙度的乘积，浪费算力。
• 逻辑错误：自身到自身的最短路径长度为0，原代码的比值计算无意义，不符合迂曲度的物理定义（实际路径长度与直线距离的比值）。

2. 连通性计算的逻辑错误 #

• 原代码将3D网格重塑为2D矩阵后计算连通分量，仅实现了每层内部的2D连通，完全忽略了上下层之间的3D连通，不符合储层3D连通性的需求。

优化后的实现代码 #

修正连通性计算（3D空间连通）+ 大幅优化迂曲度计算 #

import numpy as np
from scipy.sparse import csr_matrix, lil_matrix
from scipy.sparse.csgraph import connected_components, dijkstra

connectivity_threshold = 0.125
tortuosity_threshold = 1

def calculate_connectivity(porosity, threshold=connectivity_threshold):
    n_layers, n_rows, n_columns = porosity.shape
    total_nodes = n_layers * n_rows * n_columns
    
    # 生成有效孔隙度的掩码
    mask = porosity > threshold
    mask_flat = mask.flatten()
    
    # 构建3D邻接矩阵（LIL格式方便构建，最后转CSR）
    adj_matrix = lil_matrix((total_nodes, total_nodes), dtype=bool)
    
    # 同一层内的左右相邻
    valid_left = mask[:, :, :-1] & mask[:, :, 1:]
    left_indices = np.where(valid_left.flatten())[0]
    adj_matrix[left_indices, left_indices + 1] = True
    adj_matrix[left_indices + 1, left_indices] = True
    
    # 同一层内的上下相邻
    valid_up = mask[:, :-1, :] & mask[:, 1:, :]
    up_indices = np.where(valid_up.flatten())[0]
    adj_matrix[up_indices, up_indices + n_columns] = True
    adj_matrix[up_indices + n_columns, up_indices] = True
    
    # 层间垂直相邻
    valid_layer = mask[:-1, :, :] & mask[1:, :, :]
    layer_indices = np.where(valid_layer.flatten())[0]
    adj_matrix[layer_indices, layer_indices + n_rows * n_columns] = True
    adj_matrix[layer_indices + n_rows * n_columns, layer_indices] = True
    
    # 转CSR格式用于连通分量计算
    adj_matrix = adj_matrix.tocsr()
    
    # 计算3D连通分量
    n_components, labels = connected_components(adj_matrix, directed=False)
    
    # 还原为3D格式，非有效网格标记为0
    connectivity = np.zeros_like(porosity, dtype=int)
    connectivity[mask] = labels[mask_flat] + 1  # 组件编号从1开始
    
    return connectivity

def calculate_tortuosity(permeability_field, porosity_field, threshold=tortuosity_threshold):
    n_layers, n_rows, n_columns = permeability_field.shape
    total_nodes = n_layers * n_rows * n_columns
    
    # 预计算有效渗透率掩码
    effective_perm = permeability_field * porosity_field
    mask = effective_perm > threshold
    mask_flat = mask.flatten()
    valid_indices = np.where(mask_flat)[0]
    
    if len(valid_indices) == 0:
        return np.zeros_like(permeability_field, dtype=float)
    
    # 建立全局索引与有效索引的映射
    global_to_valid = {g_idx: v_idx for v_idx, g_idx in enumerate(valid_indices)}
    valid_to_global = {v_idx: g_idx for v_idx, g_idx in enumerate(valid_indices)}
    n_valid = len(valid_indices)
    
    # 构建有效网格的邻接矩阵
    adj_matrix = lil_matrix((n_valid, n_valid), dtype=int)
    
    # 处理同一层内的左右相邻
    valid_left = mask[:, :, :-1] & mask[:, :, 1:]
    left_global = np.where(valid_left.flatten())[0]
    for g_idx in left_global:
        v1 = global_to_valid[g_idx]
        v2 = global_to_valid[g_idx + 1]
        adj_matrix[v1, v2] = 1
        adj_matrix[v2, v1] = 1
    
    # 处理同一层内的上下相邻
    valid_up = mask[:, :-1, :] & mask[:, 1:, :]
    up_global = np.where(valid_up.flatten())[0]
    for g_idx in up_global:
        v1 = global_to_valid[g_idx]
        v2 = global_to_valid[g_idx + n_columns]
        adj_matrix[v1, v2] = 1
        adj_matrix[v2, v1] = 1
    
    # 处理层间垂直相邻
    valid_layer = mask[:-1, :, :] & mask[1:, :, :]
    layer_global = np.where(valid_layer.flatten())[0]
    for g_idx in layer_global:
        v1 = global_to_valid[g_idx]
        v2 = global_to_valid[g_idx + n_rows * n_columns]
        adj_matrix[v1, v2] = 1
        adj_matrix[v2, v1] = 1
    
    adj_matrix = adj_matrix.tocsr()
    
    # 按连通分量拆分计算
    n_components, comp_labels = connected_components(adj_matrix, directed=False)
    tortuosity_values = np.zeros_like(permeability_field, dtype=float)
    
    # 预计算所有有效节点的3D坐标
    coords = np.array([(g_idx // (n_rows * n_columns), 
                        (g_idx % (n_rows * n_columns)) // n_columns, 
                        g_idx % n_columns) 
                       for g_idx in valid_indices])
    
    for comp in range(n_components):
        comp_valid_indices = np.where(comp_labels == comp)[0]
        comp_global_indices = [valid_to_global[v_idx] for v_idx in comp_valid_indices]
        comp_size = len(comp_valid_indices)
        
        if comp_size == 1:
            # 单个节点，迂曲度设为1
            g_idx = comp_global_indices[0]
            k, i, j = g_idx // (n_rows * n_columns), (g_idx % (n_rows * n_columns)) // n_columns, g_idx % n_columns
            tortuosity_values[k, i, j] = 1.0
            continue
        
        # 提取当前分量的子邻接矩阵
        sub_adj = adj_matrix[comp_valid_indices][:, comp_valid_indices]
        
        # 计算分量内所有节点对的最短路径长度
        dist_matrix = dijkstra(sub_adj, directed=False, unweighted=True)
        
        # 找到每个节点的最大最短路径长度及对应最远节点
        max_dist = np.max(dist_matrix, axis=1)
        farthest_idx = np.argmax(dist_matrix, axis=1)
        
        # 计算欧氏直线距离
        node_coords = coords[comp_valid_indices]
        farthest_coords = node_coords[farthest_idx]
        euclid_dist = np.linalg.norm(node_coords - farthest_coords, axis=1)
        
        # 计算迂曲度（避免除零）
        tortuosity = np.where(euclid_dist > 1e-6, max_dist / euclid_dist, 1.0)
        
        # 赋值回3D数组
        for v_idx, g_idx in enumerate(comp_global_indices):
            k = g_idx // (n_rows * n_columns)
            i = (g_idx % (n_rows * n_columns)) // n_columns
            j = g_idx % n_columns
            tortuosity_values[k, i, j] = tortuosity[v_idx]
    
    return tortuosity_values

def main():
    # 生成测试数据
    porosity = np.random.rand(51, 125, 125)
    permeability = np.random.rand(51, 125, 125)
    
    # 计算连通性和迂曲度
    connectivity = calculate_connectivity(porosity, threshold=connectivity_threshold)
    tortuosity = calculate_tortuosity(permeability, porosity, threshold=tortuosity_threshold)
    
    # 保存结果
    np.save("tortuosity.npy", tortuosity)
    np.save("connectivity.npy", connectivity)

if __name__ == "__main__":
    main()

优化关键点说明 #

1. 连通性计算修正：

   • 构建真正的3D邻接矩阵，包含层内水平邻接和层间垂直邻接，实现符合需求的3D空间连通分量分析。
   • 使用稀疏矩阵存储邻接关系，大幅减少内存占用（避免79万节点稠密矩阵的600GB内存需求）。
   • 用numpy向量化操作替代三重循环，提升邻接矩阵构建效率。

2. 迂曲度计算优化：

   • 仅处理有效网格，减少计算规模。
   • 按连通分量拆分计算，避免全局处理大型图，降低单次计算复杂度。
   • 使用scipy底层C实现的dijkstra替代NetworkX的Python实现，性能提升几个数量级。
   • 修正迂曲度计算逻辑：采用最短路径长度与欧氏直线距离的比值，符合储层迂曲度的物理定义，避免原代码的除零错误。

3. 内存优化：

   • 仅保留有效节点的坐标和邻接关系，减少不必要的内存占用。
   • 使用LIL→CSR的稀疏矩阵转换，兼顾构建效率和计算效率。

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