基于元胞自动机的深部采矿岩体变形的演化模型

D01:10.13374j.isml00103x.2007.11.004 第29卷第11期 北京科技大学学报 Vol.29 No.11 2007年11月 Journal of University of Science and Technology Beijing Now.2007 基于元胞自动机的深部采矿岩体变形的演化模型 周文略杨鹏蔡嗣经 北京科技大学金属矿床高效开采与安全教有部重点实验室，北京100⑧3 摘要根据深部采矿岩体变形的离散性特征，运用元胞自动机建模方法构建深部采岩体的变形演化模型：利用元胞岩石 的状态变化概率和元胞状态变化的影响系数.实现深部采矿岩体的动力学模拟.在演化模型的基础上，进行了模拟实验，且运 用R/S分析法对结果进行分析，得到了以离散度作为岩体稳定性的标准. 关键词深部岩体变形：演化规则：元胞自动机：R/S分析；离散度 分类号TD211 近些年迅速发展起来的复杂性科学为岩体变形 岩体整体的变形破坏，通过对岩体微元状态变化分 破坏机理的研究提供了新的思路和活力，目前有关 析，建立合适的演化规则，模拟不同岩块状态的空间 研究涉及的理论和方法主要有元胞自动机(cellular 分布的演化，提供一种预测岩体变形破坏的新方法. automata简记作CA)、混沌(chaos)、非线性动力学 元胞自动机由元胞、元胞的状态空间、邻居及局 (non--linear dy namics)、自组织(sef-organization) 部规则四部分组成，可用一个四元组表示到： 等.美国数学家Hurd和Culik等人在20世纪90 A=(Ld,S,N.F). 年代初，对元胞自动机进行了严格地描述和定义：元 其中，A为一个元胞自动机：Ld为元胞空间，d为 胞自动机是一类时间、空间、格点状态都离散的动力 元胞空间的维数；S为元胞的有限状态集；N表示 系统，是一种在随机初始条件下，通过构造数学规 一个所有邻域内元胞的组合；F表示将Sm映射到S 则，来描述离散动力系统内部单元之间的自组织演 上的一个状态转换函数. 化过程的一种数学模型.使用元胞自动机方法，通 根据深部矿床岩体变形、破坏过程的岩体不同 过对构成系统的各局部之间的相互作用的建模，构 工作状态和变形的不一致条件，以及岩体结构强度 建出合适的演化规则，可以有效地模拟岩体由平衡 的不均匀性和各向异性的基本特点，本文提出了二 状态到非平衡的演化过程.这种基于元胞自动机原 维元胞自动机模型仿真的基本构造如下， 理的深部岩体的破坏机理研究，是自学习、非线性动 (1)元胞.元胞又可称为单元，元胞分布在离 态处理、演化识别、分布式表达等综合集成研究 散的二维空间的晶格点上.将巷道顶部围岩中的水 模式. 平状岩体作为研究对象，将被研究对象划分成N× 1 基于元胞自动机的岩体变形系统的 N的岩石块体每一岩石块体作为元胞.其位置用 基本模型 (x,y)表示.每个元胞的属性可以完全相同，也可 以不同 深部岩体成分各异，各种岩块强度不同，可视为 (2)元胞的状态空间.设元胞状态为元胞岩块 具有不同尺寸等级岩块构造成的块系集合.深部岩 体的工作状态集合S-{i,j}(t),表示第i行j列元 体结构的变形具有非协调、非连续特点，是一个具有 胞在时刻1的工作状态.岩块元胞体工作状态空间 复杂结构和行为特征的开放系统和耗散系统：其显 (S)包括元胞体受到外力作用但没有破裂的处于抵 著特征是高度的非线性和复杂性：简单的局部规则 抗状态（十1）、元胞体处于自然平衡状态(0)和元胞 的变化能产生整体的复杂行为，这是典型的复杂性 体处于破裂状态（一1），即为： 特征.岩体中部分岩块体的变形破坏可能会引起 S={-1,0,+1} (1) 收稿日期：2006-06-22修回日期：2006-0813 (3)邻居形式.定义二维元胞自动机的邻居分 基金项目：国家自然科学基金资助项目(N0.50374005):教有部博士 布形式为Mooe型，见图1.以最常用的规则四方 点基金资助项目(No.20020008001) 网格划分，黑色元胞为中心元胞灰色元胞为其邻 作者简介：周文略(1978一).男博士研究生：杨鹏(1965一).男. 教授博士生导师 居，以其当前状态计算中心元胞在下一时刻的状态

基于元胞自动机的深部采矿岩体变形的演化模型 周文略 杨 鹏 蔡嗣经 北京科技大学金属矿床高效开采与安全教育部重点实验室, 北京 100083 摘 要 根据深部采矿岩体变形的离散性特征, 运用元胞自动机建模方法构建深部采矿岩体的变形演化模型;利用元胞岩石 的状态变化概率和元胞状态变化的影响系数, 实现深部采矿岩体的动力学模拟.在演化模型的基础上, 进行了模拟实验, 且运 用 R/ S 分析法对结果进行分析, 得到了以离散度作为岩体稳定性的标准. 关键词 深部岩体变形;演化规则;元胞自动机;R/ S 分析;离散度 分类号 TD211 收稿日期:2006-06-22 修回日期:2006-08-13 基金项目:国家自然科学基金资助项目(No .50374005);教育部博士 点基金资助项目(No .20020008001) 作者简介:周文略(1978—), 男, 博士研究生;杨 鹏(1965—), 男, 教授, 博士生导师 近些年迅速发展起来的复杂性科学为岩体变形 破坏机理的研究提供了新的思路和活力 .目前有关 研究涉及的理论和方法主要有元胞自动机(cellular automata , 简记作CA)、混沌(chaos)、非线性动力学 (non-linear dy namics)、自 组 织 (self-organization) 等 [ 1] .美国数学家 Hurd 和 Culik 等人在 20 世纪 90 年代初,对元胞自动机进行了严格地描述和定义:元 胞自动机是一类时间 、空间、格点状态都离散的动力 系统, 是一种在随机初始条件下 , 通过构造数学规 则,来描述离散动力系统内部单元之间的自组织演 化过程的一种数学模型.使用元胞自动机方法 , 通 过对构成系统的各局部之间的相互作用的建模 , 构 建出合适的演化规则, 可以有效地模拟岩体由平衡 状态到非平衡的演化过程 .这种基于元胞自动机原 理的深部岩体的破坏机理研究 ,是自学习、非线性动 态处理 、演化识别 、分布式表达等综合集成研究 模式 . 1 基于元胞自动机的岩体变形系统的 基本模型 深部岩体成分各异, 各种岩块强度不同, 可视为 具有不同尺寸等级岩块构造成的块系集合.深部岩 体结构的变形具有非协调 、非连续特点 ,是一个具有 复杂结构和行为特征的开放系统和耗散系统;其显 著特征是高度的非线性和复杂性 ;简单的局部规则 的变化能产生整体的复杂行为 , 这是典型的复杂性 特征[ 2] .岩体中部分岩块体的变形破坏可能会引起 岩体整体的变形破坏 , 通过对岩体微元状态变化分 析 ,建立合适的演化规则,模拟不同岩块状态的空间 分布的演化 ,提供一种预测岩体变形破坏的新方法. 元胞自动机由元胞 、元胞的状态空间 、邻居及局 部规则四部分组成 ,可用一个四元组表示[ 3] : A =(Ld , S , N , F). 其中 , A 为一个元胞自动机;Ld 为元胞空间, d 为 元胞空间的维数;S 为元胞的有限状态集;N 表示 一个所有邻域内元胞的组合 ;F 表示将Sn 映射到 S 上的一个状态转换函数 . 根据深部矿床岩体变形、破坏过程的岩体不同 工作状态和变形的不一致条件, 以及岩体结构强度 的不均匀性和各向异性的基本特点, 本文提出了二 维元胞自动机模型仿真的基本构造如下. (1)元胞 .元胞又可称为单元, 元胞分布在离 散的二维空间的晶格点上.将巷道顶部围岩中的水 平状岩体作为研究对象 ,将被研究对象划分成 N × N 的岩石块体, 每一岩石块体作为元胞, 其位置用 (x , y)表示.每个元胞的属性可以完全相同, 也可 以不同 . (2)元胞的状态空间.设元胞状态为元胞岩块 体的工作状态集合 S -{i , j}(t),表示第 i 行j 列元 胞在时刻 t 的工作状态 .岩块元胞体工作状态空间 (S )包括元胞体受到外力作用但没有破裂的处于抵 抗状态(+1)、元胞体处于自然平衡状态(0)和元胞 体处于破裂状态(-1),即为 : S ={-1 , 0 , +1} (1) (3)邻居形式.定义二维元胞自动机的邻居分 布形式为 Moore 型 , 见图 1 .以最常用的规则四方 网格划分 , 黑色元胞为中心元胞, 灰色元胞为其邻 居 ,以其当前状态计算中心元胞在下一时刻的状态. 第 29 卷 第 11 期 2007 年 11 月 北 京 科 技 大 学 学 报 Journal of University of Science and Technology Beijing Vol.29 No.11 Nov.2007 DOI :10.13374/j .issn1001 -053x.2007.11.004

。1070 北京科技大学学报 第29卷 P∈(0.5,1]时，表示元胞体的强度随着围压的增大 而增大，元胞体的强度随着概率趋向于1而增大：而 当P=0.5时，表示元胞的强度无变化.具体算法如 表1所示.表1中十1表示当前元胞下一时刻的可 能状态，s为当前元胞的邻居行为，p为当前元胞 图1元胞邻居分布形式 的状态转换概率 Fig 1 Form of neighbor 表1元胞自动机演化规则 (4)规则.一般情况下，元胞下一时刻的状态 Table I Evdlutive rule of cellular automata 受其邻居元胞状态、自身状态和控制变量的影 s+1 响4，规则就是一个状态转移函数.用公式表示为： 抵抗 平衡 破裂 s1=F(品：R) 抵抗 (1-p)/2 (1-p/2 SrL=(SrL()SrL(m)) 平衡 (1-p)/2 p (1-p)/2 (2) s∈S 破裂 (1-p)/2 (1-p/2 p 1=0,1,2,… 其中，S为元胞状态空间，包括三种岩石工作状态： 实质上，表1可以概括为一个3×3的元胞状态 ,1表示元胞空间中的位置为r的元胞在十1时 转换概率矩阵，邻居元胞的状态决定当前元胞的下 一时刻的可能状态，而强度变化概率则影响当前元 刻的状态，s表示元胞空间中位置为r的元胞在t 胞的状态转换概率（转为邻居元胞的状态的概率）. 时刻的状态，sL表示位置为r的元胞的邻居L在t (2)岩体的宏观能量平衡是在一定的约束条件 时刻的状态.R是控制变量，F表示元胞自动机的 下维持的.定义这种宏观因素是以一个全局效应系 演化规则，n是邻居元胞的个数.状态变化如图2. 数来增加或减少元胞的破裂概率，设为全局系数，用 --1s-1-1+1 s广1 s壮1， s+1 M表示，取值范围为一1,1]，用系数M对原有的 -1 乐什1 》 s 强度变化概率矩阵进行调整， s+h+1,+1 广引+ 当被研究岩体整体强度高，不易出现局部变形， 图2元胞状态变化规则 则元胞体抵抗围压的概率将增加，并有M>O,在 Fig.2 Rule of cellular updated M的取值上限，强度变化概率矩阵可以表示为： P= 在上述一般规则框架下，本算法假定围岩岩体 pI1+M(p12+pB)p12(1-M)pI3(1-M) 中元胞下一个时刻的状态只受其邻居元胞的状态、 p21十M(p22+pz)p22(1-M)p2(1-M) 元胞的极限强度和全局效应因素的影响，邻居元胞 Lp31十M(p2十ps)p2(1-M)p2(1-M 不包括其自身.在这三种因素中，邻居元胞的工作 当岩体整体强度低，易出现局部变形，则元胞体 方式是状态变量函数，表现为当前时刻何种工作方 式在它们中占多数.而元胞的极限强度和全局效应 破裂的概率将增加，并有O,在M的取值下限， 因素为控制变量. 强度变化概率矩阵可以表示为： P= 为了进一步的描述元胞体的演化规则，引进下 列参数变量， pI1(1+M)p2(1+M)pB-M(p1+p12 (1)设每一个元胞体都有一个强度变化概率， P21(1+M)p2(1+M)p3-M(p21+p2) 记为P-{i,j},表示空间中第i行j列元胞的强度 Lp31(I+M)p2(1+M)P8-M(p31+p2 变化概率.岩体强度与组成岩体的岩石力学性质有 当全局无影响时，M=0,元胞体的强度变化概 关，岩石的强度会受围压、水影响以及工程施工等因 率不受影响， 素的影响发生变化，进而影响整个岩体的稳定性. 2演化模型的实现与结果分析 强度变化概率表示当前元胞以此概率转为邻居元胞 中当前时刻占多数的行为状态.P的取值范围为 2.1输入参数 [0,当P∈[0,0.5)时，元胞体的强度与围压是成 (1)网格空间大小（用G表示）.此网格空间为 反比的，元胞体强度随着概率趋向于零而降低；当 局部空间棋盘式网格，包括长和宽两个参数

图 1 元胞邻居分布形式 Fig.1 Form of neighbor (4)规则.一般情况下, 元胞下一时刻的状态 受其邻居元胞状态、自身状态和控制变量的影 响[ 4] ,规则就是一个状态转移函数 .用公式表示为 : s t +1 r =F(s t r , s t rL ;R) s t rL =(s r rL(1), …, s t rL(n)) s ∈ S t =0 , 1 , 2 , … (2) 其中, S 为元胞状态空间 , 包括三种岩石工作状态 : s t +1 r 表示元胞空间中的位置为 r 的元胞在 t +1 时 刻的状态, s t r 表示元胞空间中位置为 r 的元胞在 t 时刻的状态 , s t rL表示位置为 r 的元胞的邻居 L 在 t 时刻的状态 .R 是控制变量 , F 表示元胞自动机的 演化规则 , n 是邻居元胞的个数 .状态变化如图 2 . s t i -1, j -1 s t i -1 , j s t i -1 , j +1 s t i , j -1 s t i , j s t i , j +1 s t i +1, j -1 s t i +1 , j s t i +1 , j +1 F s t +1 i -1, j-1 s t +1 i -1, j s t+1 i -1 , j +1 s t +1 i, j -1 s t +1 i , j s t +1 i , j +1 s t +1 i +1, j-1 s t +1 i +1, j s t+1 i +1 , j +1 图 2 元胞状态变化规则 Fig.2 Rule of cellular updated 在上述一般规则框架下 , 本算法假定围岩岩体 中元胞下一个时刻的状态只受其邻居元胞的状态 、 元胞的极限强度和全局效应因素的影响, 邻居元胞 不包括其自身.在这三种因素中, 邻居元胞的工作 方式是状态变量函数, 表现为当前时刻何种工作方 式在它们中占多数.而元胞的极限强度和全局效应 因素为控制变量 . 为了进一步的描述元胞体的演化规则, 引进下 列参数变量. (1)设每一个元胞体都有一个强度变化概率 , 记为 P -{i , j},表示空间中第 i 行j 列元胞的强度 变化概率 .岩体强度与组成岩体的岩石力学性质有 关,岩石的强度会受围压 、水影响以及工程施工等因 素的影响发生变化 , 进而影响整个岩体的稳定性 . 强度变化概率表示当前元胞以此概率转为邻居元胞 中当前时刻占多数的行为状态 .P 的取值范围为 [ 0 , 1] ,当 P ∈[ 0 , 0.5)时 ,元胞体的强度与围压是成 反比的 ,元胞体强度随着概率趋向于零而降低 ;当 P ∈(0.5 , 1] 时 ,表示元胞体的强度随着围压的增大 而增大 ,元胞体的强度随着概率趋向于 1 而增大 ;而 当 P =0.5 时,表示元胞的强度无变化.具体算法如 表 1 所示.表 1 中 s t+1表示当前元胞下一时刻的可 能状态, s t L 为当前元胞的邻居行为, p 为当前元胞 的状态转换概率. 表1 元胞自动机演化规则 Table 1 Evolutive rule of cellular automata s t L s t +1 抵抗 平衡 破裂 抵抗 p (1-p)/ 2 (1-p)/ 2 平衡 (1-p)/ 2 p (1-p)/ 2 破裂 (1-p)/ 2 (1-p)/ 2 p 实质上 ,表 1 可以概括为一个 3 ×3 的元胞状态 转换概率矩阵 ,邻居元胞的状态决定当前元胞的下 一时刻的可能状态, 而强度变化概率则影响当前元 胞的状态转换概率(转为邻居元胞的状态的概率). (2)岩体的宏观能量平衡是在一定的约束条件 下维持的.定义这种宏观因素是以一个全局效应系 数来增加或减少元胞的破裂概率 ,设为全局系数 ,用 M 表示 ,取值范围为[ -1 , 1] .用系数 M 对原有的 强度变化概率矩阵进行调整 . 当被研究岩体整体强度高, 不易出现局部变形, 则元胞体抵抗围压的概率将增加 , 并有 M >0 , 在 M 的取值上限, 强度变化概率矩阵可以表示为 : P = p11+M(p12 +p13) p12(1 -M) p13(1-M) p21+M(p22 +p23) p22(1 -M) p23(1-M) p31+M(p32 +p33) p32(1 -M) p33(1-M) . 当岩体整体强度低 ,易出现局部变形 ,则元胞体 破裂的概率将增加 ,并有 M <0 ,在 M 的取值下限, 强度变化概率矩阵可以表示为: P = p11(1+M) p12(1 +M) p13 -M(p11+p12) p21(1+M) p22(1 +M) p23 -M(p21+p22) p31(1+M) p32(1 +M) p33 -M(p31+p32) . 当全局无影响时, M =0 ,元胞体的强度变化概 率不受影响 . 2 演化模型的实现与结果分析 2.1 输入参数 (1)网格空间大小(用 G 表示).此网格空间为 局部空间棋盘式网格, 包括长和宽两个参数. · 1070 · 北 京 科 技 大 学 学 报 第 29 卷

第11期 周文略等：基于元胞自动机的深部采矿岩体变形的演化模型 。1071。 (2)各个元胞强度变化概率（用P表示）.强度 分割型、随机型，如图3所示.其中黑色方框表示抵 变化概率是元胞的一种属性，它在整个模拟过程中 抗、灰色方块表示平衡、黑色方块表示破裂.对于 都保持不变.强度变化概率有随机型和同一型两 50X50的元胞空间网格，四种初始状态中各元胞状 种：随机型就是元胞空间中每个元胞的强度变化概 态的数量（即岩体的体积）和分配比例如下： 率都可能不同，用随机数产生每个元胞的强度变化 中心对称型：{289,800,1411}，{0.12,0.32， 概率；同一型是指元胞空间中每个元胞的强度变化 0.56: 概率相同. 横向分割型：{800,850,850，{032,0.34， (3)全局效应系数M(M∈[一1，) 0.34: (4)各个元胞的初始状态（用S:表示）.本文设 纵向分割型：{800.850,850，{032,0.34， 定四种初始状态，即中心对称型、横向分割型、纵向 0.34; 50 50 g ·抵抗 抵抗 ·平衡 平衡 40 ,抵抗 40 抵抗 30 30 20 20 10 10 20 30 40 10 203040 50 元胞空间宽度 元胞空间宽度 (a)中心对称型 ()横向分料型 抵抗 抵抗 平衡 平衡 抗 ·抵抗 30 20 10 10 20 30 40 50 10 20 30 40 50 元胞空间宽度 元胞空间宽度 (©)纵向分割型 (d)随机型 图3岩体初始状态分布的四种情形 Fig.3 Four initial arrangements of cellular rock 随机分布型：在这种初始状态中，元胞空间中的 的元胞的数量总和.元胞容量用N表示，抵抗、平 各种工作状态呈随机分布，各种状态的数量总和大 衡、破裂三种状态的容量分别为N、Nb和N. 致相等，但其具体分布和具体的数量是不确定的 (2)岩体的稳定性系数：设抵抗、平衡、破裂三 2.2输出参数 种元胞状态容量的数学期望及期望的极差分别是 输入参数可用一个四元向量来表示：{G,p, E、Eb、E:和R,如果定义岩体的离散度(D,=R/ M,S}.实验中，对每个元胞的强度变化概率在[0， (E十Eb十E),则离散度越低岩体越稳定，离散度 1]之间随机地取值，并在0，]上服从均匀分布，设 越高岩体越不稳定. {G,p,M,S}={50X50,随机，0，s},Is为元胞的 (3)岩体的最终状态：元胞状态容量变化曲线 初始状态.考查在每个元胞强度变化概率随机分布 的直观审视达到平衡，也即各个元胞态容量曲线达 时这个岩体的状态变化，即岩体的稳定性. 到总体水平时. (1)元胞状态容量：被研究岩体中在特定时刻 2.3结果分析 处于某一种状态的元胞数量的总和，如岩体中破裂 重标极差分析(R/S分析法)是赫斯特发明的

(2)各个元胞强度变化概率(用 P 表示).强度 变化概率是元胞的一种属性, 它在整个模拟过程中 都保持不变 .强度变化概率有随机型和同一型两 种:随机型就是元胞空间中每个元胞的强度变化概 率都可能不同, 用随机数产生每个元胞的强度变化 概率;同一型是指元胞空间中每个元胞的强度变化 概率相同 . (3)全局效应系数 M(M ∈[ -1 , 1] ). (4)各个元胞的初始状态(用 S i 表示).本文设 定四种初始状态 , 即中心对称型、横向分割型、纵向 分割型、随机型,如图 3 所示 .其中黑色方框表示抵 抗 、灰色方块表示平衡 、黑色方块表示破裂.对于 50 ×50 的元胞空间网格, 四种初始状态中各元胞状 态的数量(即岩体的体积)和分配比例如下 : 中心对称型:{289 , 800 , 1411}, {0.12 , 0.32 , 0.56}; 横向分割型:{800 , 850 , 850}, {0.32 , 0.34 , 0.34}; 纵向分割型:{800 , 850 , 850}, {0.32 , 0.34 , 0.34}; 图 3 岩体初始状态分布的四种情形 Fig.3 Four initial arrangements of cellular rock 随机分布型 :在这种初始状态中,元胞空间中的 各种工作状态呈随机分布, 各种状态的数量总和大 致相等,但其具体分布和具体的数量是不确定的. 2.2 输出参数 输入参数可用一个四元向量来表示:{G , p , M , S}.实验中, 对每个元胞的强度变化概率在[ 0 , 1] 之间随机地取值, 并在[ 0 , 1] 上服从均匀分布 , 设 {G , p , M , S }={50 ×50 ,随机 , 0 , Is}, Is 为元胞的 初始状态 .考查在每个元胞强度变化概率随机分布 时这个岩体的状态变化, 即岩体的稳定性. (1)元胞状态容量:被研究岩体中在特定时刻 处于某一种状态的元胞数量的总和 ,如岩体中破裂 的元胞的数量总和.元胞容量用 N 表示 , 抵抗、平 衡 、破裂三种状态的容量分别为 N r 、N b 和 N f . (2)岩体的稳定性系数:设抵抗、平衡 、破裂三 种元胞状态容量的数学期望及期望的极差分别是 E r 、E b 、Ef 和 R , 如果定义岩体的离散度(Ds =R / (E r +E b +E f)),则离散度越低岩体越稳定, 离散度 越高岩体越不稳定 . (3)岩体的最终状态 :元胞状态容量变化曲线 的直观审视达到平衡 , 也即各个元胞态容量曲线达 到总体水平时. 2.3 结果分析 重标极差分析(R /S 分析法)是赫斯特发明的. 第 11 期 周文略等:基于元胞自动机的深部采矿岩体变形的演化模型 · 1071 ·

。1072· 北京科技大学学报 第29卷 重标极差等于观测值的极差除以标准差用重标极 极差.定义： 差随观测时间的变化情况，即用赫斯特指数可以度 R (4) 量上述有偏随机游动的趋势9.重标极差(R/S)与 D,=E,十Eh+Ei 观测时间的关系如下： 表2变形系统特征值 (R/S)N=CN (3) Table 2 System eigenvalues of deformation 其中，N为观测次数或时间间隔，C为常数H为赫 抵抗状态 平衡状态 破裂状态 数组 斯特指数. H. H H 赫斯特指数在0,1)之间取值，由不同的观测次 079 129 079 405 078 14 数N可以得到一系列R/S值，以ln(N)和n(R/ 2 079 35 080 55 079 58 S)作线性回归分析，随着观测次数的增大，斜率会 3 0.80 122 080 52 080 59 逐渐增大，最终会达到一个极大值，此时对应的斜率 0.79 40 0 79 71 080 33 和观测次数即为时间序列的赫斯特指数和平均循环 0.79 35 36 080 35 长度.在重标极差分析中，涉及到平均循环长 0.78 078 17 度L. 7 0.82 30 080 50 从图4可以直接看出，关系曲线在N=29的地 8 079 34 080 357 方发生突变，此时平均循环长度L=29,对应L≤29 078 079 079 134 时H值为083. 10 0.79 54 079 88 079 91 0.84 平均值0.79 57.4079868079 848 0.83 0.82 图5为中心对称型初始状态的元胞行为岩体容 0.81 量图(T为演化步数)，其他三种与此相似直观上 -0.80 无明显区别， 0.79 0.78 950 N.....N.---N 0.77 900k 070 40 6080100120 850 800 图4赫斯特指数与观测次数关系 750 10 20 30 4050 Fig.4 Relation of Hurst index to observation number 演化步骤，T 从10组由随机初始出发的元胞状态容量，来估 图5三种元胞容量图 算随机分布强度变化概率时的赫斯特指数(H)和平 Fig.5 Capability of three types of cells 均循环长度(L).其中，每组状态容量中包括了演 对于图3所示的岩体演化空间的不同初始状 化2000步的数据.表2为10组HL的数据显示. 态，经过演化后，都从直观审视中达到平衡状态（所 表2中H的变化很小，在0.78与0.82之间：而L 谓平衡状态是指各行为的状态容量的期望值不再变 的变化很大，在14与405之间.本文实验的10组 化之时，可以在状态容量曲线上进行直观审视).计 数据，分别出自不同的强度变化概率的分布，即这 算抵抗、平衡、破裂这三种元胞状态容量的数学期望 10组数据具有不同的初始参数，因此系统特征有一 及期望的极差如下：对于中心对称初始状态，E= 定的差异.对于赫斯特指数与平均循环长度，均取 82230,Eb=820.60,E=857.10,R=3650:对于 这10组数据的平均值作为它们的估计值. 横向分割初始状态，E=835.10,Eb=831.70,E= 元胞状态容量变化曲线的赫斯特指数反映：强 833.20,R=3.40:对于纵向分割初始状态，E,= 度变化概率随机分布的岩体系统中，其岩体变形演 79840,E=843.50,E=858.10,R=59.70:对于 变遵循有偏随机游动，岩体变形是有趋势的 随机分割初始状态，E,=821.20,E=820.80,E= 定义离散度(Ds)来反映对象岩体的稳定 85800,R=37.20. 性$y.设抵抗、平衡、破裂这三种元胞行为岩体容 根据以上数学期望值和极差，由式(4)计算出离 量的数学期望分别为E、Eb、E。R为E、E、E的 散度(D,):中心对称初始状态，D=1.46%;横向分

重标极差等于观测值的极差除以标准差, 用重标极 差随观测时间的变化情况, 即用赫斯特指数可以度 量上述有偏随机游动的趋势[ 5] .重标极差(R /S )与 观测时间的关系如下 : (R / S)N =CN H (3) 其中 , N 为观测次数或时间间隔 , C 为常数, H 为赫 斯特指数 . 赫斯特指数在[ 0 , 1)之间取值 ,由不同的观测次 数 N 可以得到一系列 R/ S 值, 以 ln(N)和 ln(R/ S )作线性回归分析 , 随着观测次数的增大, 斜率会 逐渐增大 ,最终会达到一个极大值 ,此时对应的斜率 和观测次数即为时间序列的赫斯特指数和平均循环 长度 [ 6-7] .在重标极差分析中, 涉及到平均循环长 度 L . 从图 4 可以直接看出 ,关系曲线在 N =29 的地 方发生突变,此时平均循环长度 L =29 ,对应 L ≤29 时 H 值为 0.83 . 图 4 赫斯特指数与观测次数关系 Fig.4 Relation of Hurst index to observation number 从 10 组由随机初始出发的元胞状态容量 ,来估 算随机分布强度变化概率时的赫斯特指数(H)和平 均循环长度(L).其中, 每组状态容量中包括了演 化2 000 步的数据.表 2 为10 组 H-L 的数据显示 . 表2 中 H 的变化很小, 在 0.78 与 0.82 之间;而 L 的变化很大 ,在 14 与 405 之间.本文实验的 10 组 数据, 分别出自不同的强度变化概率的分布, 即这 10 组数据具有不同的初始参数 ,因此系统特征有一 定的差异.对于赫斯特指数与平均循环长度 , 均取 这 10 组数据的平均值作为它们的估计值. 元胞状态容量变化曲线的赫斯特指数反映:强 度变化概率随机分布的岩体系统中 ,其岩体变形演 变遵循有偏随机游动 ,岩体变形是有趋势的. 定义离 散度 (Ds)来反 映对 象岩体 的稳 定 性[ 8-9] .设抵抗、平衡、破裂这三种元胞行为岩体容 量的数学期望分别为 Er 、E b 、E f , R 为E r 、Eb 、Ef 的 极差.定义 : Ds = R Er +E b +E f (4) 表 2 变形系统特征值 Table 2 System eigenvalues of deformation 数组 抵抗状态 平衡状态 破裂状态 Hr L r Hb L b H f Lf 1 0.79 129 0.79 405 0.78 14 2 0.79 35 0.80 55 0.79 58 3 0.80 122 0.80 52 0.80 59 4 0.79 40 0.79 71 0.80 33 5 0.79 35 0.79 36 0.80 35 6 0.78 20 0.79 73 0.78 17 7 0.82 80 0.81 30 0.80 50 8 0.79 24 0.79 34 0.80 357 9 0.78 35 0.79 24 0.79 134 10 0.79 54 0.79 88 0.79 91 平均值 0.79 57.4 0.79 86.8 0.79 84.8 图 5 为中心对称型初始状态的元胞行为岩体容 量图(T 为演化步数), 其他三种与此相似, 直观上 无明显区别 . 图 5 三种元胞容量图 Fig.5 Capability of three types of cells 对于图 3 所示的岩体演化空间的不同初始状 态 ,经过演化后,都从直观审视中达到平衡状态(所 谓平衡状态是指各行为的状态容量的期望值不再变 化之时, 可以在状态容量曲线上进行直观审视).计 算抵抗 、平衡 、破裂这三种元胞状态容量的数学期望 及期望的极差如下 :对于中心对称初始状态 , E r = 822.30 , Eb =820.60 , E f =857.10 , R =36.50 ;对于 横向分割初始状态 , E r =835.10 , E b =831.70 , Ef = 833.20 , R =3.40 ;对于纵向分割初始状态 , E r = 798.40 , Eb =843.50 , E f =858.10 , R =59.70 ;对于 随机分割初始状态 , E r =821.20 , E b =820.80 , Ef = 858.00 , R =37.20 . 根据以上数学期望值和极差 ,由式(4)计算出离 散度(Ds):中心对称初始状态, Ds =1.46 %;横向分 · 1072 · 北 京 科 技 大 学 学 报 第 29 卷

第11期 周文略等：基于元胞自动机的深部采矿岩体变形的演化模型 。1073。 割初始状态，D=014%:纵向分割初始状态，D,= 破坏引起的地质灾害提供了一种新的方法. 239%:随机分割初始状态，D,=1.49%. 如果把10%作为离散度度量岩体稳定性的阀 参考文献 值10，那么从这四种初始状态出发的岩体都是稳定 【刂成思危。复杂性科学探索.北京：民主与建设出版社.1998 【习许强、黄润秋.非线性科学理论在地质灾害评价预测中的应 的；且离散度越低越趋向稳定. 用一地质灾害系统分析原理.山地学报.2000,18(3)：272 3结论 3 Bandini,Mauri.Serra.Celular automata from atheoretical par allel computationl model to its appication to complex systems. (1)对于复杂组成的岩体可以运用元胞自动 Parallel Comput.2001(27):539 机原理准确地模拟岩体变形演化情况。建立足够精 【4应尚军，蔡嗣经.基于元胞自动机的股票市场投资行为模拟. 确的反映深部岩体诸特点的变形、动力和运动过程 系统工程学报，2001,16(5)：382 [5 Mandelbrot BB.Statistical methodology for nom periodic cycles: 及其数学模型. from the covariance to R/S analysis.Ann Econ Social Meas (2)通过对模拟结果的重标极差分析，对岩体 1972(1):257 的变形演化进行趋势预测，并把岩体的离散度作为 【(杨正瓴，林孔元余贻鑫用时间序的雅普诺夫指数算预报电力 岩体稳定性的衡量标准，离散度越低的岩体越稳定， 系统中的某些失稳现象中国电机工程学报.2001,21(1)：5 离散度越高的岩体越不稳定.这为地质灾害的预防 [7 Wolfram S.Stat istical mecharics of cellular automata.Rev Mod Phs1983,55(3):601 提供科学的理论支持. [8 Capcarrere,Mathieu.A statistical study of a class of celular evo (3)只要对深部岩体取样，获取特定尺寸、初始 lutionary akorithms.Evdl Comput 1999,7(3):34 强度以及岩石地质分布等数据，就可以利用该基本 【9蒙立军，背树芳，王运敏.岩石断裂演化过程元胞自动机模型 模型进行演化计算，这为研究岩体的非线性动力学 的研究.金属矿山.2002(5)：58 特征提供了一种有力的数值模拟手段，为预测岩体 【10钱鸣高，石平五.矿山压力与岩层控制。北京：中国矿业大学 出版社.2003：54 Cellular automaton model for rock deformation in deep mining ZHOU Wenlue,YANG Peng,CAI Sijing Key Lab of Chinese Miristry of Educat ion for High-Efficient Mining and Safety of Metal Depsit Uriversity of Science and Techmlogy Beijing.Beijing 100083.China ABSTRACT According to discrete characteristics of deformation.an evolution model of deformation in mine rock-mass was built by using cellular automaton(CA).In order to measure the complex process of deformation in mine-rock body,the change probability of cellular rock and the influence coefficient were defined.Simulat ion tests were carried out and the experimental results were analyzed by use of R/S analysis.It is show n that dis- crete degree can be taken as the criterion of rock-mass stability. KEY WORDS deep rock-mass deformation;evolution rule;cellular automata model;R/S analy sis;discrete degree

割初始状态 , Ds =0.14 %;纵向分割初始状态 , Ds = 2.39 %;随机分割初始状态, Ds =1.49 %. 如果把 10 %作为离散度度量岩体稳定性的阀 值[ 10] , 那么从这四种初始状态出发的岩体都是稳定 的;且离散度越低,越趋向稳定 . 3 结论 (1)对于复杂组成的岩体, 可以运用元胞自动 机原理准确地模拟岩体变形演化情况, 建立足够精 确的反映深部岩体诸特点的变形 、动力和运动过程 及其数学模型. (2)通过对模拟结果的重标极差分析 , 对岩体 的变形演化进行趋势预测, 并把岩体的离散度作为 岩体稳定性的衡量标准, 离散度越低的岩体越稳定 , 离散度越高的岩体越不稳定.这为地质灾害的预防 提供科学的理论支持 . (3)只要对深部岩体取样, 获取特定尺寸、初始 强度以及岩石地质分布等数据 , 就可以利用该基本 模型进行演化计算 ,这为研究岩体的非线性动力学 特征提供了一种有力的数值模拟手段, 为预测岩体 破坏引起的地质灾害提供了一种新的方法 . 参 考 文 献 [ 1] 成思危.复杂性科学探索.北京:民主与建设出版社, 1998 [ 2] 许强, 黄润秋.非线性科学理论在地质灾害评价预测中的应 用———地质灾害系统分析原理.山地学报, 2000 , 18(3):272 [ 3] Bandini , Mau ri, S erra.Cellular automata, from a theoretical par￾allel computational model to its appli cation to complex systems. Parallel Comput , 2001(27):539 [ 4] 应尚军, 蔡嗣经.基于元胞自动机的股票市场投资行为模拟. 系统工程学报, 2001 , 16(5):382 [ 5] Mandelbrot B B .Statistical methodology f or non-peri odic cycles: from the covariance to R/ S analysis.Ann Econ Social Meas, 1972(1):257 [ 6] 杨正瓴, 林孔元, 余贻鑫.用时间序的雅普诺夫指数算预报电力 系统中的某些失稳现象.中国电机工程学报, 2001 , 21(1):5 [ 7] Wolfram S.St atistical mechani cs of cellular automata.Rev Mod Phys, 1983 , 55(3):601 [ 8] Capcarrere , Mathieu.A st atistical study of a class of cellular evo￾lutionary algorithms.Evol Comput, 1999 , 7(3):34 [ 9] 蒙立军, 肖树芳, 王运敏.岩石断裂演化过程元胞自动机模型 的研究.金属矿山, 2002(5):58 [ 10] 钱鸣高, 石平五.矿山压力与岩层控制.北京:中国矿业大学 出版社, 2003:54 Cellular automaton model for rock deformation in deep mining ZHOU Wenlue , Y ANG Peng , CAI S ijing Key Lab of Chinese Ministry of Education for High-Efficient Mining and Safety of Metal Deposit, Uni versit y of Science and Tech nology Beijing , Beijing 100083 , China ABSTRACT According to discrete characteristics of defo rmation , an evolution model of defo rmation in mine rock-mass w as built by using cellular automaton (CA).In o rder to measure the complex process of deformation in mine-rock body , the change probability of cellular rock and the influence coefficient w ere defined .Simulation tests w ere carried out and the experimental results w ere analyzed by use of R /S analysis.It is show n that dis￾crete degree can be taken as the criterion of rock-mass stability . KEY WORDS deep rock-mass deformation ;evolution rule ;cellular automata model ;R /S analy sis ;discrete deg ree 第 11 期 周文略等:基于元胞自动机的深部采矿岩体变形的演化模型 · 1073 ·