Python 自动分析岩石破坏过程:裂纹提取、扩展演化与损伤区域可视化(含代码)

引言

传统裂纹分析通常依赖人工观察破坏照片，或者通过软件手动标注裂纹位置。这种方式在裂纹数量较少时尚可满足需求，但面对复杂破裂过程时，人工统计往往存在效率低、主观性强以及难以进行连续阶段对比等问题。

随着数字图像处理技术的发展，利用 Python 可以将岩石破坏图像转化为可量化的裂纹数据。通过图像增强、阈值分割、形态学处理和骨架提取等方法，可以自动识别裂纹区域，并进一步计算裂纹长度、裂纹面积、裂纹方向以及损伤集中区域。本文以模拟岩石破坏过程图像为例，展示 Python 在裂纹自动分析中的应用流程。代码实现了从原始破坏图像输入，到裂纹提取、裂纹演化统计和损伤区域可视化的完整过程。该方法不仅适用于岩石试验照片，也可以用于 PFC/FDEM 模拟裂纹结果、混凝土破坏图像以及复合材料损伤分析。通过 Python 自动处理，可以快速建立从图像到裂纹参数的分析流程。

import osimport numpy as npimport pandas as pdimport matplotlibmatplotlib.use("Agg")import matplotlib.pyplot as pltfrom scipy.ndimage import gaussian_filterfrom skimage import (    morphology,    filters,    measure,    exposure,    color)np.random.seed(2026)out_dir = "crack_damage_analysis_results"os.makedirs(out_dir, exist_ok=True)plt.rcParams["font.family"] = "Times New Roman"plt.rcParams["axes.unicode_minus"] = FalseH = 600W = 900# 生成模拟岩石破坏图像def generate_crack_image(stage):    img = np.ones((H, W))*0.65    # 岩石纹理    texture = gaussian_filter(        np.random.normal(0,1,(H,W)),        sigma=15    )    texture = texture / texture.std()    img += texture*0.05    crack = np.zeros((H,W),dtype=bool)    # 主剪切裂纹    length = 180 + stage*60    x0 = 120    y0 = 500    for i in range(length):        x = x0+i*2        y = int(            y0 -            0.55*i +            20*np.sin(i/30)        )        width = 1+stage        if (            5<y<H-5 and            5<x<W-5        ):            crack[                y-width:y+width,                x-width:x+width            ] = True    # 次生裂纹    n_branch = stage*3+2    for i in range(n_branch):        start_x = np.random.randint(            200,650        )        start_y = np.random.randint(            150,450        )        angle = np.random.uniform(            -1.5,1.5        )        for j in range(            80+stage*20        ):            x = int(                start_x+j*np.cos(angle)            )            y = int(                start_y+j*np.sin(angle)            )            if (                5<x<W-5 and                5<y<H-5            ):                crack[y-1:y+2,x-1:x+2]=True    crack = morphology.binary_dilation(        crack,        morphology.disk(stage)    )    img[crack]=0.15    img += np.random.normal(        0,        0.015,        img.shape    )    img=np.clip(        img,        0,        1    )    return img, crack# 2. 裂纹自动提取函数def crack_detection(img):    # 对比增强    img_eq = exposure.equalize_adapthist(        img,        clip_limit=0.03    )    # 裂纹暗区域    threshold = filters.threshold_otsu(        img_eq    )    crack_mask = img_eq < threshold*0.92    # 去噪    crack_mask = morphology.remove_small_objects(        crack_mask,        min_size=20    )    crack_mask = morphology.binary_closing(        crack_mask,        morphology.disk(2)    )    # 骨架    skeleton = morphology.skeletonize(        crack_mask    )    return crack_mask,skeleton# 图1 裂纹自动提取img, true_crack = generate_crack_image(    stage=5)crack_mask, skeleton = crack_detection(    img)fig,ax=plt.subplots(    1,3,    figsize=(15,5))ax[0].imshow(    img,    cmap="gray")ax[0].set_title(    "Original failure image")ax[1].imshow(    crack_mask,    cmap="gray")ax[1].set_title(    "Automatic crack extraction")overlay=np.dstack(    [        img,        img,        img    ])overlay[skeleton]=[    1,    0,    0]ax[2].imshow(    overlay)ax[2].set_title(    "Crack skeleton")for a in ax:    a.axis("off")plt.tight_layout()plt.savefig(    os.path.join(        out_dir,        "fig1_crack_extraction.png"    ),    dpi=500,    bbox_inches="tight")plt.close()# 图2 裂纹扩展演化stages=np.arange(    1,    7)lengths=[]areas=[]damage=[]for s in stages:    img_s, crack_s = generate_crack_image(        s    )    sk = morphology.skeletonize(        crack_s    )    lengths.append(        sk.sum()    )    areas.append(        crack_s.sum()    )    damage.append(        crack_s.sum()/(H*W)    )df=pd.DataFrame(    {        "stage":stages,        "crack_length":lengths,        "crack_area":areas,        "damage_variable":damage    })df.to_csv(    os.path.join(        out_dir,        "crack_evolution.csv"    ),    index=False)fig,ax=plt.subplots(    1,3,    figsize=(15,4.5))ax[0].plot(    stages,    lengths,    "-o",    linewidth=2)ax[0].set_title(    "Crack length evolution")ax[0].set_xlabel(    "Loading stage")ax[0].set_ylabel(    "Skeleton pixels")ax[1].plot(    stages,    areas,    "-o",    color="darkred",    linewidth=2)ax[1].set_title(    "Crack area evolution")ax[2].plot(    stages,    damage,    "-o",    color="purple",    linewidth=2)ax[2].set_title(    "Damage variable evolution")for a in ax:    a.grid(alpha=0.2)plt.tight_layout()plt.savefig(    os.path.join(        out_dir,        "fig2_crack_evolution.png"    ),    dpi=500,    bbox_inches="tight")plt.close()# 图3 裂纹方向 + 损伤热图label_img = measure.label(    skeleton)regions = measure.regionprops(    label_img)angles=[]for r in regions:    if r.area<10:        continue    angles.append(        np.degrees(            r.orientation        )    )angles=np.array(    angles)fig,ax=plt.subplots(    1,2,    figsize=(12,5))# --------方向玫瑰图--------theta=np.deg2rad(    angles)bins=np.linspace(    0,    np.pi,    19)hist,_=np.histogram(    theta,    bins)centers=(bins[:-1]+bins[1:])/2ax0=plt.subplot(    121,    polar=True)ax0.bar(    centers,    hist,    width=np.pi/18,    alpha=0.75)ax0.set_title(    "Crack orientation distribution")# --------裂纹密度热图--------density=gaussian_filter(    crack_mask.astype(float),    sigma=25)ax1=plt.subplot(    122)im=ax1.imshow(    density,    cmap="hot")ax1.set_title(    "Crack density map")ax1.axis("off")plt.colorbar(    im,    ax=ax1,    shrink=0.8)plt.tight_layout()plt.savefig(    os.path.join(        out_dir,        "fig3_crack_direction_damage.png"    ),    dpi=500,    bbox_inches="tight")plt.close()

特别声明：

以上代码与文案均为网上资料整合而成，仅供广大同行们参考学习，如有侵权请联系删除。