638 工程科学学报，第43卷，第5期 同的矿物胶结而成，不同矿物若以均质材料进行 中长石、石英和云母所占图像面积百分比，分别 赋值校准很可能导致不准确的模拟结果.可以基 为49.85%、29.74%和20.41%，其中云母含量过高， 于真实岩石材料数字图像、以矿物颗粒表面不同 这与阈值选取有关.实际花岗岩矿物组分中除上 颜色或亮度变化准确识别花岗岩表面细观结构， 述三种主要矿物外，还有角闪石、少量辉石等次要 进而将所得结构作为建立离散元模型的基本条件 矿物，可将长石和石英以外的矿物百分比视为黑 试验岩块尺寸为：50mm×50mm×100m,使 云母矿物含量 用数码相机获得花岗岩试块的数字图像，图像中 的每个像素点都可以用像素点矩阵表示： Gray value 0 f(1,1)f1,2) f1,N) Feldspar f2,1) f2,2) f2,N) F(i,)= (1) Quartz Mica fM,1)f(M,2) f(M,N) 255 式中，f亿，)表示在0,255)范围内的像素点(，) 的灰度值，M和N代表像素点的行数和列数.采 图1花岗岩图像矿物识别.(a)标准试件：(b)局部放大图 Fig.I Mineral recognition from a granite image:(a)standard test 用MATLAB中图像处理方法，对原始图像进行平 specimen;(b)partial enlarged detail 滑处理以消除噪声影响.然后采用阈值分割算法， 采用图像处理技术得到的花岗岩细观组分， 识别花岗岩矿物组分特征. 并不能完全展现其表面形态特征.在二维数值计 改进的最大类间方差法适用于图像的多阈值 算模型中，通常以圆形颗粒为基础来构建矿物组 分割，以类间特性差最大或类内特性差最小为判 分的分布区域，从而建立反映试件组分分布的细 据，求得图像分割的最佳阈值.设灰度范围O, 观模型：然而，矿物组分的三维形态特征难以直接 L]的图像中，第i级像素点个数为N,像素点总个 获取，PFC3D中的球形颗粒很难用于直接建立表征 数为N=方N,则第1级像素点的出现概率为 矿物组分分布的多面体.为简化计算模型，本文根 =0 据试样二维图像，通过编程识别不同矿物组分，计 N 选择灰度阈值组T={,2,…,tm-1小，将图像 算矿物组分占比，进而建立相应的三维分析模型. 分割为m个区域，灰度区域的灰度值分别为 1.2细观参数标定 C1=[0,t,C2=[t1,2],,Cm=[tm-1,L,则方差最 在岩石微观结构研究中21-2四，通常认为矿物成 大值的最优阈值T为： 分的强度控制着裂纹的萌生，长石、石英、黑云母 (T")=arg ma ,- 矿物强度关系为1：1.17：0.42，压痕试验下长石、 (2) 石英、黑云母的断裂韧度比值为1：2.08：0.76.根 据不同矿物破坏应力阈值的裂纹演化特征，发现 =∑Pd,=∑iP:lwx:X= 3) 微裂纹在黑云母中首先产生、在石英中产生最晚) icC. ieC, 在岩石宏细观参数关系研究中以长石组分颗 其中，o(s=1,2,3,…,m)为每一灰度区域像素点占 粒的黏结强度为平行黏结力学参数基本值，石英 图像总比的比例，C,为区域灰度值的集合，入为平 和云母组分颗粒的黏结强度取基本值的1.6倍和 均灰度值，为图像的总体灰度.采用该算法进行 0.3倍，不同组分的黏结强度取相应组分的平均 双阈值求解，得到灰度图像的最佳分割阈值为 值，颗粒的最小半径Rmn取1.2mm,最大半径和最 124和168，从而对花岗岩中的长石、石英、云母图 小半径比为1.66 像进行分割.本文尚未考虑图像中不同矿物组分 模型细观力学性质参数需要根据室内试验结 之间存在灰度交叉的情况，由此导致的误差予以 果进行标定，以保证模型试样与实验试样的宏观 忽略 力学响应基本一致.本文使用AW2000型全数字 花岗岩标准试样的灰度图如图1(a)所示，岩 电液控制刚性压力试验机进行完整岩石单轴压缩 石结构特征和阈值分割如图1(b)所示.通过对灰 试验，位移控制加载速率为0.03 mm:min、加载至 度图像中各像素点进行数量统计，得到矿物成分 试样破坏.室内单轴压缩试验得到的宏观力学性同的矿物胶结而成，不同矿物若以均质材料进行 赋值校准很可能导致不准确的模拟结果. 可以基 于真实岩石材料数字图像、以矿物颗粒表面不同 颜色或亮度变化准确识别花岗岩表面细观结构， 进而将所得结构作为建立离散元模型的基本条件. 试验岩块尺寸为：50 mm × 50 mm × 100 m，使 用数码相机获得花岗岩试块的数字图像，图像中 的每个像素点都可以用像素点矩阵表示： F(i, j) =   f(1,1) f(1,2) ··· f(1,N) f(2,1) f(2,2) ··· f(2,N) . . . . . . . . . f(M,1) f(M,2) ··· f(M,N)   （1） 式中，f (i，j) 表示在 [0, 255] 范围内的像素点 (i，j) 的灰度值，M 和 N 代表像素点的行数和列数. 采 用 MATLAB 中图像处理方法，对原始图像进行平 滑处理以消除噪声影响. 然后采用阈值分割算法， 识别花岗岩矿物组分特征. N = ∑ L i=0 Ni Pi = Ni N T = {t1,t2,··· ,tm−1} C1 = [0,t1] C2 = [t1,t2], ··· Cm = [tm−1,L] δ 2 B T ∗ 改进的最大类间方差法适用于图像的多阈值 分割，以类间特性差最大或类内特性差最小为判 据，求得图像分割的最佳阈值. 设灰度范围 [0， L] 的图像中，第 i 级像素点个数为 Ni，像素点总个 数 为 ， 则 第 i 级像素点的出现概率为 . 选择灰度阈值组 ，将图像 分 割 为 m 个区域 ，灰度区域的灰度值分别为 ， ， ，则方差 最 大值的最优阈值 为： δ 2 B (T ∗ ) = arg max   ∑m s=1 ωs(λs −λ) 2   （2） ωs = ∑ i∈Cs Pi ;λs = ∑ i∈Cs iPi/ωs ;λ = ∑ L i=0 iPi （3） ωs(s = 1,2,3,··· ,m) Cs λs λ 其中， 为每一灰度区域像素点占 图像总比的比例， 为区域灰度值的集合， 为平 均灰度值， 为图像的总体灰度. 采用该算法进行 双阈值求解 ，得到灰度图像的最佳分割阈值为 124 和 168，从而对花岗岩中的长石、石英、云母图 像进行分割. 本文尚未考虑图像中不同矿物组分 之间存在灰度交叉的情况，由此导致的误差予以 忽略. 花岗岩标准试样的灰度图如图 1（a）所示，岩 石结构特征和阈值分割如图 1（b）所示. 通过对灰 度图像中各像素点进行数量统计，得到矿物成分 中长石、石英和云母所占图像面积百分比，分别 为 49.85%、29.74% 和 20.41%，其中云母含量过高， 这与阈值选取有关. 实际花岗岩矿物组分中除上 述三种主要矿物外，还有角闪石、少量辉石等次要 矿物，可将长石和石英以外的矿物百分比视为黑 云母矿物含量. (a) (b) Feldspar Mica Quartz 255 0 Gray value 图 1    花岗岩图像矿物识别. （a）标准试件；（b）局部放大图 Fig.1     Mineral  recognition  from  a  granite  image:  (a)  standard  test specimen; (b) partial enlarged detail 采用图像处理技术得到的花岗岩细观组分， 并不能完全展现其表面形态特征. 在二维数值计 算模型中，通常以圆形颗粒为基础来构建矿物组 分的分布区域，从而建立反映试件组分分布的细 观模型；然而，矿物组分的三维形态特征难以直接 获取，PFC3D 中的球形颗粒很难用于直接建立表征 矿物组分分布的多面体. 为简化计算模型，本文根 据试样二维图像，通过编程识别不同矿物组分，计 算矿物组分占比，进而建立相应的三维分析模型. 1.2    细观参数标定 在岩石微观结构研究中[21−22] ，通常认为矿物成 分的强度控制着裂纹的萌生，长石、石英、黑云母 矿物强度关系为 1∶1.17∶0.42，压痕试验下长石、 石英、黑云母的断裂韧度比值为 1∶2.08∶0.76. 根 据不同矿物破坏应力阈值的裂纹演化特征，发现 微裂纹在黑云母中首先产生、在石英中产生最晚[23] . 在岩石宏细观参数关系研究中[24] ，以长石组分颗 粒的黏结强度为平行黏结力学参数基本值，石英 和云母组分颗粒的黏结强度取基本值的 1.6 倍和 0.3 倍，不同组分的黏结强度取相应组分的平均 值，颗粒的最小半径 Rmin 取 1.2 mm，最大半径和最 小半径比为 1.66. 模型细观力学性质参数需要根据室内试验结 果进行标定，以保证模型试样与实验试样的宏观 力学响应基本一致. 本文使用 AW2000 型全数字 电液控制刚性压力试验机进行完整岩石单轴压缩 试验，位移控制加载速率为 0.03 mm·min−1、加载至 试样破坏. 室内单轴压缩试验得到的宏观力学性 · 638 · 工程科学学报，第 43 卷，第 5 期