书 书 书 第５１卷 第２期 ２０１６年４月 西 南 交 通 大 学 学 报 ＪＯＵＲＮＡＬ ＯＦ ＳＯＵＴＨＷＥＳＴ ＪＩＡＯＴＯＮＧ ＵＮＩＶＥＲＳＩＴＹ Ｖｏｌ． ５１ Ｎｏ． ２ Ａｐｒ． ２０１６ 收稿日期：２０１５１２０１ 基金项目：国家自然科学基金资助项目（４１１７２３２１，４１５７１００４）；国家重点基础研究发展计划课题（２００８ＣＢ４２５８０２） 引文格式：姚令侃，黄艺丹． 山地系统灾变行为自组织临界性研究［Ｊ］． 西南交通大学学报，２０１６，５１（２）：３１３３３０． 文章编号：０２５８２７２４（２０１６）０２０３１３１８ ＤＯＩ：１０． ３９６９ ／ ｊ． ｉｓｓｎ． ０２５８２７２４． ２０１６． ０２． ０１１
山地系统灾变行为自组织临界性研究
姚
令
侃
１，２，３，
黄
艺
丹
１，２，３ （１． 西南交通大学土木工程学院，四川成都６１００３１；２． 高速铁路线路工程教育部重点实验室，四川成都 ６１００３１；３． 抗震工程技术四川省重点实验室道路与铁道工程抗震技术研究所，四川成都６１００３１） 摘 要：为了研究山地系统宏观动力学的整体规律，基于戴维斯地貌发育理论，提出处于地貌发育阶段幼年晚 期的河谷及壮年期的山地具有自组织临界性（ｓｅｌｆ ｏｒｇａｎｉｚｅｄ ｃｒｉｔｉｃａｌｉｔｙ，ＳＯＣ）的内禀属性，并建立了基于斯特拉勒 积分的大流域地貌发育阶段判别方法．以地震触发崩塌滑坡为切入点，通过震区实震资料分析、元胞自动机模 拟、振动台沙堆模型实验，提出不同烈度区地震触发崩塌滑坡分布的演化规律：在Ⅶ度、Ⅷ度、Ⅸ度区，崩塌滑坡 规模与出现频率之间存在良好的负幂律关系，在Ⅹ度区，幂律关系弱化，在Ⅺ度区，这一关系转为对数正态分布． 通过３个案例介绍了ＳＯＣ理论在岩土体地震扰动深度评估、泥石流防治工程设计径流量极值计算、基于地震活 动性参数ｂ值在地应力评估中的应用． 关键词：山地系统；自组织临界性；地震触发崩塌滑坡；元胞自动机 中图分类号：Ｐ６４２ 文献标志码：Ａ ＳｅｌｆＯｒｇａｎｉｚｅｄ Ｃｒｉｔｉｃａｌｉｔｙ ｏｆ Ｍｏｕｎｔａｉｎ Ｓｙｓｔｅｍ Ｃａｔａｓｔｒｏｐｈｉｃ Ｂｅｈａｖｉｏｒｓ ＹＡＯ Ｌｉｎｇｋａｎ１，２，３， ＨＵＡＮＧ Ｙｉｄａｎ１，２，３ （１． Ｓｃｈｏｏｌ ｏｆ Ｃｉｖｉｌ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｓｏｕｔｈｗｅｓｔ Ｊｉａｏｔｏｎｇ Ｕｎｉｖｅｒｓｉｔｙ，Ｃｈｅｎｇｄｕ ６１００３１，Ｃｈｉｎａ；２． ＭＯＥ Ｋｅｙ Ｌａｂｏｒａｔｏｒｙ ｏｆ ＨｉｇｈＳｐｅｅｄ Ｒａｉｌｗａｙ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｃｈｅｎｇｄｕ ６１００３１，Ｃｈｉｎａ；３． Ｒｏａｄ ａｎｄ Ｒａｉｌｗａｙ Ｅｎｇｉｎｅｅｒｉｎｇ Ｒｅｓｅａｒｃｈ Ｉｎｓｔｉｔｕｔｅ， Ｓｉｃｈｕａｎ Ｋｅｙ Ｌａｂｏｒａｔｏｒｙ ｏｆ Ｓｅｉｓｍｉｃ Ｅｎｇｉｎｅｅｒｉｎｇ ａｎｄ Ｔｅｃｈｎｏｌｏｇｙ，Ｃｈｅｎｇｄｕ ６１００３１，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｉｎ ｏｒｄｅｒ ｔｏ ｓｔｕｄｙ ｔｈｅ ｍａｃｒｏｄｙｎａｍｉｃ ｂｅｈａｖｉｏｒ ｏｆ ｍｏｕｎｔａｉｎ ｓｙｓｔｅｍｓ，ａ ｖｉｅｗ ｐｏｉｎｔ ｗａｓ ｐｒｏｐｏｓｅｄ ｂｙ Ｄａｖｉｓｓ ｔｈｅｏｒｙ ｏｆ ｔｈｅ ｅｒｏｓｉｏｎ ｃｙｃｌｅ ｔｈａｔ ｖａｌｌｅｙｓ ｉｎ ｌａｔｅ ｃｈｉｌｄｈｏｏｄ ｏｒ ｍｏｕｎｔａｉｎｓ ｉｎ ｍｉｄｄｌｅ ａｄｕｌｔｈｏｏｄ ｈａｖｅ ｉｎｔｒｉｎｓｉｃ ｐｒｏｐｅｒｔｉｅｓ ｏｆ ｓｅｌｆｏｒｇａｎｉｚｅｄ ｃｒｉｔｉｃａｌｉｔｙ（ＳＯＣ）． Ｂａｓｅｄ ｏｎ Ｓｔｒａｈｌｅｒｓ ｉｎｔｅｇｒａｌ，ａ ｍｅｔｈｏｄ ｆｏｒ ｄｉｓｔｉｎｇｕｉｓｈｉｎｇ ｌａｎｄｆｏｒｍ ｄｅｖｅｌｏｐｍｅｎｔ ｓｔａｇｅｓ ｉｎ ａ ｌａｒｇｅ ｂａｓｉｎ ｉｓ ｅｓｔａｂｌｉｓｈｅｄ． Ｆｏｒ ｌａｎｄｓｌｉｄｅｓ ｔｒｉｇｇｅｒｅｄ ｂｙ ｅａｒｔｈｑｕａｋｅ，ｔｈｅ ｄａｔａ ｏｆ ｌａｎｄｓｌｉｄｅｓ ｃｏｌｌｅｃｔｅｄ ｉｎ ｅａｒｔｈｑｕａｋｅ ｚｏｎｅｓ ｗｅｒｅ ａｎａｌｙｚｅｄ，ａ ｃｅｌｌｕｌａｒ ａｕｔｏｍａｔａ ｍｏｄｅｌ ｗａｓ ｂｕｉｌｔ，ａｎｄ ｓｈａｋｉｎｇ ｔａｂｌｅ ｅｘｐｅｒｉｍｅｎｔｓ ｗｅｒｅ ｐｅｒｆｏｒｍｅｄ． Ａｓ ａ ｒｅｓｕｌｔ，ａｎ ｅｖｏｌｕｔｉｏｎ ｌａｗ ｏｆ ｌａｎｄｓｌｉｄｅｓ ｄｉｓｔｒｉｂｕｔｉｏｎ ｉｎ ｄｉｆｆｅｒｅｎｔ ｓｅｉｓｍｉｃ ｉｎｔｅｎｓｉｔｙ ｚｏｎｅｓ ｗａｓ ｒｅｖｅａｌｅｄ． Ｓｐｅｃｉｆｉｃａｌｌｙ，ｔｈｅｒｅ ｅｘｉｓｔｓ ａ ｓｔｒｏｎｇ ｎｅｇａｔｉｖｅ ｐｏｗｅｒｌａｗ ｒｅｌａｔｉｏｎｓｈｉｐ ｂｅｔｗｅｅｎ ｔｈｅ ｓｉｚｅｓ ａｎｄ ｆｒｅｑｕｅｎｃｉｅｓ ｏｆ ｌａｎｄｓｌｉｄｅｓ ｉｎ ｚｏｎｅｓ ｗｉｔｈ ｓｅｉｓｍｉｃ ｉｎｔｅｎｓｉｔｙ Ⅶ，Ⅷ ａｎｄ Ⅸ；ｔｈｅ ｒｅｌａｔｉｏｎｓｈｉｐ ｂｅｃｏｍｅｓ ａ ｗｅａｋ ｐｏｗｅｒ ｌａｗ ｉｎ ｚｏｎｅｓ ｗｉｔｈ ｓｅｉｓｍｉｃ ｉｎｔｅｎｓｉｔｙ Ⅹ，ａｎｄ ｃｈａｎｇｅｓ ｉｎｔｏ ａ ｌｏｇｎｏｒｍａｌ ｄｉｓｔｒｉｂｕｔｉｏｎ ｉｎ ｚｏｎｅｓ ｗｉｔｈ ｓｅｉｓｍｉｃ ｉｎｔｅｎｓｉｔｙ Ⅺ． Ｆｉｎａｌｌｙ， ａｐｐｌｉｃａｔｉｏｎｓ ｏｆ ＳＯＣ ｔｏ ｓｏｉｌ ｄｅｐｔｈ ａｓｓｅｓｓｍｅｎｔ ｕｎｄｅｒ ｅａｒｔｈｑｕａｋｅ ｄｉｓｔｕｒｂａｎｃｅ，ｅｘｔｒｅｍｅ ｒｕｎｏｆｆ ｃａｌｃｕｌａｔｉｏｎ ｉｎ ｄｅｂｒｉｓ ｆｌｏｗ ｃｏｎｔｒｏｌ ｄｅｓｉｇｎ，ａｎｄ ｃｒｕｓｔａｌ ｓｔｒｅｓｓ ｅｖａｌｕａｔｉｏｎ ｂａｓｅｄ ｏｎ ｔｈｅ ｓｅｉｓｍｉｃ ａｃｔｉｖｉｔｙ ｐａｒａｍｅｔｅｒ ｂｖａｌｕｅ，ａｒｅ ｉｎｔｒｏｄｕｃｅｄ ｔｈｒｏｕｇｈ ｔｈｒｅｅ ｃａｓｅ ｓｔｕｄｉｅｓ． Ｋｅｙ ｗｏｒｄｓ：ｍｏｕｎｔａｉｎ ｓｙｓｔｅｍ；ｓｅｌｆｏｒｇａｎｉｚｅｄ ｃｒｉｔｉｃａｌｉｔｙ（ＳＯＣ）；ｌａｎｄｓｌｉｄｅｓ ｔｒｉｇｇｅｒｅｄ ｂｙ ｅａｒｔｈｑｕａｋｅ；西 南 交 通 大 学 学 报 第５１卷 ｃｅｌｌｕｌａｒ ａｕｔｏｍａｔａ 在非线性物理学中，作为灾变理论提出的前沿 理论———自组织临界状态（ｓｅｌｆｏｒｇａｎｉｚｅｄ ｃｒｉｔｉｃａｌｉｔｙ， ＳＯＣ）理论，是Ｐ Ｂａｋ等为解释无序的、非线性复杂 系统的行为特征而提出的新概念．这类系统包含着 众多发生短程相互作用的组元，并自发地向着一种 临界状态进化．在临界状态下，小事件引起的连锁 反应能对系统中大量数目的组元产生影响，从而导 致大规模事件的发生．虽然发生的小事件比大事件 多，但遍及所有规模的连锁反应是动态特性的一个 必不可少的部分，所有时空关联函数都是幂律 （ｐｏｗｅｒｌａｗ）的，因此幂律可以作为自组织临界状 态的证据［１］ ． 沙堆模型是ＳＯＣ 的范例． Ｇ Ａ Ｈｅｌｄ等采用在 圆盘上逐粒加沙的方式构造沙堆，当沙堆的倾角在 临界角附近时沙堆停止增长，此时，对新添加沙粒 的响应是无法预测的，沙粒可能固定在沙堆上，也 可能引起小范围沙粒的滑动，还可能导致更大规模 的雪崩（ａｖａｌａｎｃｈｅ），但总是呈现崩塌规模与出现 频率成反比的幂律关系［２］ ． ＳＯＣ的概念已被用于解释从山脉形成到股市 波动等许多复杂现象．地震学是ＳＯＣ 最早关注的 领域，如古登堡里克特定律表明地震的频度与震 级之间存在幂律关系． Ｐ Ｂａｋ认为这种幂律来源于 地壳系统的ＳＯＣ，反之，也可把幂律作为地壳被锁 定在永久临界状态的证据［３］ ． 在地学领域研究 ＳＯＣ是具有学科交叉性质的前沿研究，从宏观层 面，於崇文院士以完整和独立的命题提出了固体地 球系统的复杂性与自组织临界性［４］，认为地质系 统是自然界中的一种异常复杂的开放、远离平衡、 相互作用的巨大耗散动力系统，具有自组织临界性 的内禀属性，其时空行为服从地质作用的自组织临 界过程动力学．从微观层面，许强等［５］论述了岩石 破裂过程的自组织临界特征． 运用ＳＯＣ理论认识和描述山地系统灾变动力 学行为是本团队的主攻方向．国家自然科学基金委 先后立项资助本团队开展了“散粒体自组织临界 性”（２００３）、“大尺度散粒体自组织临界性及其判 据”（２００５）、“地震触发崩塌滑坡自组织临界性” （２０１２）项目研究，这些题目也反映了本团队从散 粒体材料层面至山地系统层面ＳＯＣ性质的研究历 程．本文主要以地震触发崩塌滑坡作为研究山地系 统灾变行为ＳＯＣ的切入点，开展一项专题研究，并 介绍在应用领域开展的一些探索． １
山
地
系
统
具
备
ＳＯＣ性
质
的
条
件
１． １ ＳＯＣ系统的必要条件 从现象学的角度，ＳＯＣ系统应具有的必要条件 是［６］： （１）系统是耗散的，包含大量（数目应在数百 万个以上）发生短程相互作用的组元；大量自由度 以某种均衡态势存在，不出现优势自由度而使系统 仅有几个集体自由度来表征的现象． （２）组元之间最近邻位置相互作用，但存在着 长程相关的关系． （３）系统自发地朝着临界状态进化，并将会永 久性地或在一个有意义的时段内被锁定在这个 状态． （４）系统同时具备两种动力学特征：一是敏感 性，系统处于临界状态，微小的扰动也可能产生遍 及全局的连锁反应；二是鲁棒性，系统的整体不受 破坏． （５）所有的时空关联函数均为幂次． 从必要条件可看出，ＳＯＣ系统的基本特征是以 临界状态为动力学吸引子． 即在外界持续的物质 （或能量）供给条件下，广延耗散动力系统自发地 向临界状态演化，一旦到达临界状态，将通过与外 界能量的交换，永久性地或在一个有意义的时段内 能被锁定在这个状态．这是从自然系统演化趋势角 度提出的ＳＯＣ 基本条件，其特点是强调自然系统 ＳＯＣ的时段性． 但以上５条尚不能涵盖ＳＯＣ 系统需要满足的 所有条件，还需要补充一些能够引起系统产生自组 织到临界态的趋势的特定特征． 从构造性的角度，ＳＯＣ 行为只会在缓慢驱动 的、相互作用占主导地位的阈值系统中出现．阈的 存在使系统积累能量、应力和物质直至达到不稳定 阈值． ＳＯＣ 是慢驱动机制，即能量积累速率缓慢， 而能量耗散的时间尺度远小于前者．以地震为例， 地壳处于运动中，因而经受缓慢的形变，形变将导 致岩石中应力的积累．如果在某些地方岩石不能再 承受应力，就会发生突发事件，并向周围环境释放 能量．如果环境可以承担这种能量的增加，那么这 是一个孤立事件；如果环境也处在一种临界状态， 就可能发生链式反应，地震就相当于应力释放的链 ４ １ ３第２期 姚令侃，等：山地系统灾变行为自组织临界性研究 式反应，涉及到比较大的区域．这两种过程的时间 标度差别巨大，如“５·１２”汶川大地震，专家估计应 力在龙门山推覆构造带上的积累过程可能长达千 年，地震时应力释放的主要过程仅在约１００ ｓ内完 成［７］ ． １． ２ 基于戴维斯地貌发育理论的山地系统临界状 态思辨 美国著名地理学家戴维斯（Ｗ Ｍ Ｄａｖｉｓ）于 １８９９ 年首次创立了侵蚀循环学说（ｔｈｅｏｒｙ ｏｆ ｔｈｅ ｃｙｃｌｅ ｏｆ ｅｒｏｓｉｏｎ），认为地块从开始上升到被逐渐剥 蚀夷平，直至降低到起伏不大的地面或者接近基准 面的准平原之间，存在着连续的、同时又有阶段性 的剥蚀过程和地表形态．在地表发育的过程中， Ｄａｖｉｓ强调构造、作用和时间（侵蚀阶段）这３个要 素之间的相互作用影响［８］，进而将循环过程中的 地形发展分为：地形起伏不大，河间地广阔平坦的 幼年期；地面主要由谷坡和狭窄的分水岭组成的壮 年期；具有残丘的准平原的老年期（图１）． （ａ）最初，地形起伏和 缓，流水不畅 狭（ｂ窄）幼，高年地早宽期阔，沟平坦缘 主（ｃ，仍）幼有年沟晚缘期，平，岩坦坡高为地 （ｄ）壮年期，多为岩坡 与狭窄的分水岭 起（ｅ伏）壮较缓年晚，谷期底，地宽形展 （ｆ蚀）老余残年山期的，成准为平具原有 （第ｇ二）再循次环，重构造现抬幼升年，进早入期 图１ 侵蚀循环示意［９］ Ｆｉｇ． １ Ｓｃｈｅｍａｔｉｃ ｏｆ ｅｒｏｓｉｏｎ ｃｙｃｌｅ［９］ 按照侵蚀循环理论，地貌发育循环过程中的地 形变化如图２所示，图中，曲线Ｌ１、Ｌ２ 分别代表原 始地面河谷的平均高度和分水岭的平均高度随时 间的变化．虚线表示地壳短暂快速的抬升阶段，Ｈ１ 表示原始地形的起伏程度． 由图２可知，整个幼年期阶段直至壮年中期， 河谷强烈迅速下切，河谷平均海拔持续下降，河谷 之间由比较宽广的分水地发展为分水岭，平均海拔 高度没有显著变化，因此，原始地面的起伏程度迅 速增加，至壮年中期，起伏程度达到最大，分水岭与 河谷之间的高差达到最大值Ｈ２．壮年期晚期和老 年期阶段，河谷已经接近侵蚀基准面，高程基本没 有变化，分水岭受到外营力的侵蚀而不断降低，分 水岭变得浑圆低矮，直至准平原．
图２ Ｄａｖｉｓ河流地貌发育图式（根据Ｗ． Ｍ． Ｄａｖｉｓ） Ｆｉｇ． ２ Ｔｈｅ ｓｅｑｕｅｎｃｅ ｉｎ ｔｈｅ ｄｅｖｅｌｏｐｍｅｎｔａｌ ｃｈａｎｇｅｓ ｏｆ ｌａｎｄｆｏｒｍｓ ｏｆ ａｎ ｉｄｅａｌ ｇｅｏｇｒａｐｈｉｃａｌ ｃｙｃｌｅ ｂｙ Ｗ． Ｍ． Ｄａｖｉｓｓ ｔｈｅｏｒｙ 地貌的形成和发展是内、外营力相互作用的结 果．内营力趋向于使山体隆升，增强区域起伏程度， 即使山体愈发陡峻；外营力趋向于使山体高度降 低、削平，减弱区域起伏程度，即使山体愈发浑圆低 ５ １ ３
西 南 交 通 大 学 学 报 第５１卷 矮．内营力是系统的，来自全球性的板块构造，可以 认为其速度是确定的，而外营力是变化的．在新构 造运动强烈的山区，山地隆升速度每年可达几毫米 至十几毫米，山体剥蚀降低的速度小于地块隆升速 度；但随着山体高度的上升，带来侵蚀基准面的降 低，进而加大地表物质的重力作用和水流的下切侵 蚀与搬运作用，外营力效应随之增强．外营力效应 绝对量的增加必将导致内外营力作用的相对差距 缩小，当地貌演化处于内外营力相当的阶段时，由 于二者对山地地貌塑造的反向效应，山体坡度最大 只能达到一个特定值，即所谓的临界坡度（对应于 图２的Ｈ２ 时期），按照ＳＯＣ理论，临界坡度成为该 阶段山地斜坡系统演化的吸引子． 对照图１，虽然在一个侵蚀循环内，地貌将经 历幼年期、壮年期、老年期，流域内的坡体相应经历 向临界坡度发展、达到临界坡、偏离临界坡的演变 过程，但处于地貌发育阶段幼年晚期地块的河谷 （图１（ｃ）），以及壮年期的山地（图１（ｄ）），其斜坡 系统能够保持在临界坡度，即在这一时段内已经演 化到了临界状态． 若从构造性的角度，在漫长的地质过程中，地 表形态的变化是组成地表的物质在渐变与突变的 不断转化中进行的．渐变是指地形在漫长时间完成 的变化过程，突变则是在急促的短时间内完成的变 化过程．由于板块运动，在长达千百万年的时间尺 度上使山脉隆升，将能量（势能）缓慢而持续地输 入斜坡系统，使斜坡系统能量积累到临界状态，这 一过程是非随机的；而在外营力随机扰动条件下， 斜坡物质可能失稳并在重力作用下向下输移，即斜 坡系统能量又通过突发性的山地灾害的方式耗散， 这一过程是随机的．山地隆升驱动斜坡系统构造的 速度（ｃｍ ／ ａ）和崩塌滑坡滑动速度（ｍ ／ ｓ）之间出现 时间标度的巨大分离，相互作用为主导、阈动力学 和慢驱动的有机结合，使得斜坡系统的演化成为自 然界中ＳＯＣ的良好实例． 综上所述，无论从现象学的角度，还是从构造 性的角度，对照ＳＯＣ的基本定义，可以认为当且仅 当处于地貌发育阶段幼年晚期的河谷，如横断山三 江并流区的高山峡谷地段（图３、图４），以及壮年 期的山地，如龙门山（图５），山地斜坡系统成为这 些区域地貌的主体，斜坡系统已经演化到了临界状 态，ＳＯＣ成为系统的内禀属性，它的灾变行为将服 从ＳＯＣ过程动力学．特别是崩塌、滑坡、泥石流等 斜坡重力作用类地表过程，其共同特征是能量耗散 以斜坡物质失稳下滑实现的，从物理概化角度，可 认为沙堆模型是一种能抓住真实斜坡系统基本特 征的简单理想模型，因而它们的动力学特征都应能 在ＳＯＣ 的概念框架下得到解释． 图３ 三江并流区怒江河谷已演化到临界状态 Ｆｉｇ． ３ Ｎｕｊｉａｎｇ ｒｉｖｅｒ ｖａｌｌｅｙｓ ｉｎ ｔｈｅ ｔｈｒｅｅ ｐａｒａｌｌｅｌ ｒｉｖｅｒｓ ｒｅｇｉｏｎ ｈａｖｉｎｇ ｅｖｏｌｖｅｄ ｔｏ ｔｈｅ ｃｒｉｔｉｃａｌ ｓｔａｔｅ 图４ 三江并流区澜沧江河谷已演化到临界状态 Ｆｉｇ． ４ Ｌａｎｃａｎｇ ｒｉｖｅｒ ｖａｌｌｅｙｓ ｉｎ ｔｈｒｅｅ ｐａｒａｌｌｅｌ ｒｉｖｅｒｓ ｒｅｇｉｏｎ ｈａｖｉｎｇ ｅｖｏｌｖｅｄ ｔｏ ｔｈｅ ｃｒｉｔｉｃａｌ ｓｔａｔｅ 图５ 龙门山已演化到临界状态， 在强震作用下大面积崩塌 Ｆｉｇ． ５ Ｌｏｎｇｍｅｎ ｍｏｕｎｔａｉｎｓ ｈａｖｉｎｇ ｅｖｏｌｖｅｄ ｔｏ ｔｈｅ ｃｒｉｔｉｃａｌ ｓｔａｔｅ，ｗｉｔｈ ｌｏｔｓ ｏｆ ｌａｎｄｓｌｉｄｅｓ ｔｒｉｇｇｅｒｅｄ ｂｙ ｅａｒｔｈｑｕａｋｅ １． ３ 基于斯特拉勒积分的山地系统发育阶段判别 方法 以上提出了以地貌发育阶段作为界定山地系 统具有ＳＯＣ性质的适用范围，应用中还必须建立 地貌发育阶段的定量判别方法． １９５２年，斯特拉勒（Ａ Ｎ Ｓｔｒａｈｌｅｒ）提出侵蚀流 ６ １ ３
第２期 姚令侃，等：山地系统灾变行为自组织临界性研究 域的面积高程分析方法［１０］，可以定量推求 Ｄａｖｉｓ 的地貌发育阶段．面积高程曲线分析法是描述一 定高度范围内的面积随相对高度变化所表示的曲 线及其所围成的面积，相对高度可以确定侵蚀过程 的强度，而残留的面积可以代表这种强度下地貌的 保持能力，因此，可以说面积高程曲线提供了地貌 的发育信息［１１］ ． 记流域内等高线的值和最低点之间的高差为 ｈ，每条等高线以上的面积为ａ，全流域面积为Ａ０， 流域内最高点和最低点之间的高差为Ｈ，分别以 Ｘ ＝ ａ ／ Ａ０， （１） 为横坐标和纵坐标画图，可以得到曲线 Ｙ ＝ ｈ ／ Ｈ， （２） 也就是面积高程曲线，称为斯特拉勒曲线（ｔｈｅ Ｓｔｒａｈｌｅｒｓ ｃｕｒｖｅ）． 设定积分 Ｓ ＝ ∫ １． ０ ０ ｆ（ｘ）， （３） 式中： Ｓ为斯特拉勒曲线与坐标轴包围的面积，称为 斯特拉勒积分（ｔｈｅ Ｓｔｒａｈｌｅｒｓ ｉｎｔｅｇｒａｌ）． 用Ｓ值推求侵蚀流域地貌演化阶段，即： Ｓ ＞ ０． ６０，幼年期； ０． ３５≤Ｓ≤０． ６０，壮年期； Ｓ ＜ ０． ３５，老年期． 由于斯特拉勒积分为无量纲参数，因此，曲线 可以描述和比较不同规模的流域，但必须是在流域 内地貌处于同一发育阶段的前提下，所以，斯特拉 勒积分应用于小流域一般没有问题．然而，对于大 型山区河流，不同河段可能处于不同的地貌单元， 这与Ｄａｖｉｓ循环理论中流域内的各个部分同时发 育的假定是相悖的，不能直接应用该方法． 作者提出包含不同地貌单元的大流域斯特拉 勒积分计算方法： 首先确定计算区域，通过选择河道控制点（也 就是计算流域的出口点），确定适用于Ｄａｖｉｓ 侵蚀 循环理论的范围．选择河道控制点的原则为： ① 控制点以上的流域是灾势评估的研究区； ② 河道控制点以上的流域应基本上处于同一 地貌单元． 然后利用ＤＥＭ和ＡｒｃＧＩＳ技术，获得斯特拉勒 曲线的横、纵坐标值．具体方法是在具有规则格网 的ＤＥＭ上，对研究区域内所有大于某个高程值的 栅格单元个数进行累加，再乘以栅格单元的面积Ａ０， 即得到该高程值以上的面积［１２］，如式（ ４）所示． Ａｈ０ ＝ ∑ ｈｍａｘ ｈ ＝ ｈ０ ＮｈＡ０， （４） 式中： Ａｈ０为等高线值为ｈ０ 以上的面积； Ｎｈ为研究区域内高程值为ｈ的栅格个数． 以（Ａｈ０／ Ａ，ｈ０／ Ｈ）为横、纵坐标值，绘制斯特拉 勒曲线，可计算得到斯特拉勒积分． 现举例说明基于斯特拉勒积分的山地系统发 育阶段判别方法． ２０１０年４月１４ 日，我国青海省 玉树县发生Ｍｓ７． １ 级地震，影响流域主要水系包 括通天河、扎曲、巴曲等，地形以高海拔、低起伏为 主．震区位于通天河和金沙江的交汇地段，圈定的 研究区域的流域包括沱沱河、通天河、金沙江上游 在内的长江上游河段，研究区域的水系格局如图６ 所示． 图６ 研究流域及周边流域、 裂点位置、计算流域区域图 Ｆｉｇ． ６ Ｒｅｇｉｏｎ ｇｒａｐｈｓ ｏｆ ｔｈｅ ｓｔｕｄｙ ｂａｓｉｎ ａｎｄ ｉｔｓ ｓｕｒｒｏｕｎｄｉｎｇ ｂａｓｉｎ，ｔｈｅ ｋｎｉｃｋ ｐｏｉｎｔ ｌｏｃａｔｉｏｎ， ａｎｄ ｔｈｅ ｃｏｍｐｕｔａｔｉｏｎ ｂａｓｉｎ 由于青藏高原的分阶段、非均一隆升，使高原 东缘的外流河产生溯源侵蚀，并形成新的河谷，它 与未被溯源侵蚀的老河谷交替的地方，河床坡度突 然增加，形成裂点［１３］，裂点上下游的河谷往往处于 不同的地貌发育状态．通过提取河道的纵剖面图 （图７）可以看出，在海拔４ ０００ ｍ 的地方存在裂 点．目前青藏高原普遍存在着高原面、盆地面两级 夷平面．盆地面是青藏高原形成于上新世纪初至上 新世纪末的第二级夷平面，盆地面内平坦开阔，切 割微弱．位于裂点以上的金沙江上游谷地是溯源侵 蚀尚未达到的地方，这一地区的盆地面保存完整， 较少切割．裂点以下的金沙江河谷位于我国地形最 为陡峻的横断山区，河流切割成深邃的峡谷．裂点 ７ １ ３
西 南 交 通 大 学 学 报 第５１卷 上、下游两个地貌单元差异明显［１４］ ．又因玉树地区 处于该裂点的上方，可满足上文提出的选择河道控 制点的两条原则，故决定选择该裂点为计算流域的 出口点，绘制斯特拉勒曲线，计算斯特拉勒积分 Ｓ ＝ ０． ３２，如图８所示．图８ 表明，震区内流域属于 壮年晚期发育阶段，宽谷缓丘是主要地貌特征．在 玉树地震时，地震触发的崩塌滑坡无论从规模还是 数量上，都远小于汶川地震和芦山地震，认为地貌 发育阶段不同是玉树地震崩塌滑坡灾害小于汶川、 芦山地震的主要因素（汶川和芦山地震区均处于 壮年期）．
图７ 研究河流纵剖面图 Ｆｉｇ． ７ Ｌｏｎｇｉｔｕｄｉｎａｌ ｐｒｏｆｉｌｅ ｏｆ ｔｈｅ ｓｔｕｄｙ ｒｉｖｅｒ 图８ 计算区域的斯特拉勒曲线 Ｆｉｇ． ８ Ｓｔｒａｈｌｅｒ ｃｕｒｖｅ ｏｆ ｔｈｅ ｃｏｍｐｕｔａｔｉｏｎ ｂａｓｉｎ 反之，若按图６所示的全河段计算相对应流域 的斯特拉勒积分，得Ｓ′ ＝ ０． ５３，据此推断整个流域 的发育期处于壮年期，应该是坡陡谷深的地貌，显 然与震区的自然景观不符．至此，解决了斯特拉勒 分析方法只能适应于小流域的问题，建立了山地系 统发育阶段的定量判别方法． ２不同烈度区地震触发崩塌滑坡分
布
规
律
发生在四川龙门山地区的“５·１２”汶川地震、 “４·２０”芦山地震，是人类有现代观测仪器以来，地 震触发崩塌滑坡数量最多、资料最翔实的两次大地 震，震后有关地震触发崩塌滑坡的研究成为热点课 题，但国内外大多数研究集中在对该场地震崩塌滑 坡灾害的统计、地震触发崩塌滑坡形成的力学机理 等方面，缺少从物理学层面对其总体特征进行描述 的整体理论． 龙门山是青藏高原边缘山脉中的陡度变化最 大的山脉，在３０多公里范围内海拔从７００ ｍ升高 到５ ０００ ｍ．晚新生代中新世以来，龙门山至少有 ５ ～ １０ ｋｍ 的底层被剥蚀掉，上升速度约达 ０． ６ ｍｍ ／ ａ．近年来的地形变化资料表明，该构造带 的九顶山地区正以０． ３ ～ ０． ４ ｍｍ ／ ａ速度持续上 升［１５］，这种隆升和夷平的持续作用，造成龙门山河 谷深切、地势陡峻的地貌景观．龙门山整体处于地 貌演化的壮年早期，斜坡系统已经演化到了临界状 态，这就是能在ＳＯＣ 的概念框架下，研究汶川地 震、芦山地震触发崩塌滑坡整体分布规律的前提条 件．为此，选择地震触发崩塌滑坡作为研究山地系 统灾变行为ＳＯＣ的切入点，在原型问题的物理代 表性、资料拥有量等方面，均具有独特的优势． ２． １ 地震触发崩塌滑坡实震资料分析 利用震后遥感影像资料进行人工目视解译是 大面积获取地震触发崩塌滑坡信息的主要方法，由 于使用的遥感影像资料精度不同、判识人员的判识 标准和经验不同等，对同一区域的判识可能会出现 较大差异，因此，现场调查工作不可忽视．但以上细 节并不是关键问题，影响地震触发崩塌滑坡分布规 律最具控制性的因素是地震烈度，在研究崩塌滑坡 分布规律时，应该按地震烈度区分别统计． 汶川地震后，获取２００８ 年６ 月４ 日的ＡＬＯＳ 卫星影像（分辨率为１０ ｍ），覆盖范围约为 ９ ７５０ ｋｍ２（ ３０°５８． ７８′Ｎ ～ ３２°２． ３′Ｎ，１０３°３３． ９７′Ｅ ～ １０４°３６． １７′Ｅ），包括北川、汶川、茂县、都江堰等区 域．另外，还通过下载获得汶川震前ＴＭ卫星图片． 由于这套ＡＬＯＳ遥感数据资料Ⅸ度以下烈度区云 层覆盖率高，因此，仅对Ⅸ度及以上烈度区进行崩 塌滑坡的遥感解译工作．在上述区域内，采用人机 交互的目视解译技术，辅以野外实地调查，判译出 地震触发崩塌滑坡９ ３４１ 处，解译结果如 图９所示． 芦山地震后，通过多种途径收集和购买了震区 航空、航天遥感影像资料，主要包括中国科学院遥 感与数字地球研究所提供的三批航片图（拍摄时 间为２０１３年４月２０ 日和２１日，其中，第一批影像 ８ １ ３第２期 姚令侃，等：山地系统灾变行为自组织临界性研究 分辨率为０． ６ ｍ，第二、三批影像有０． ４ ｍ 和 ２ ｍ两种分辨率）；购买的资源三号卫星影像（拍摄 时间为２０１３年５月１３ 日，分辨率为２． １ ｍ）．对以 上遥感影像资料进行几何纠正、融合、拼接、图像增 强处理后，获得芦山震区遥感数据５ ６５５ ｋｍ２ （１０２°２７． ６８′Ｅ ～ １０３° ３０． ８４′ Ｅ，２９° ４４． ５４′ Ｎ ～ ３０° ２６． ９１′Ｎ），包括邛崃市、名山县、雨城区、芦山县、 天全县、宝兴县、荥经县等区域．在上述区域，对地 震触发的崩塌滑坡灾害点进行逐一解译，并结合实 地考察检验解译结果，最终确定地震触发崩塌滑坡 １ ６０８处，解译结果如图１０所示． 图９ 汶川地震触发的崩塌滑坡分布示意图 Ｆｉｇ． ９ Ｔｈｅ ｄｉｓｔｒｉｂｕｔｉｏｎ ｍａｐ ｏｆ ｌａｎｄｓｌｉｄｅｓ ｉｎｄｕｃｅｄ ｂｙ Ｗｅｎｃｈｕａｎ ｅａｒｔｈｑｕａｋｅ 按不同烈度区，分别对汶川Ⅸ度、Ⅹ度、Ⅺ度， 芦山Ⅶ度、Ⅷ度、Ⅸ度区的崩塌滑坡进行统计分析， 以崩塌滑坡面积（Ａ）作为规模的度量，分析崩塌滑 坡规模与发生频率之间的关系．由于芦山地震各烈 度区的遥感影像资料更为完备（Ⅸ度区遥感覆盖 面积为１００％，Ⅷ度区遥感覆盖面积为９４． ５％，Ⅶ 度区遥感覆盖面积为５９． ４％），故芦山震区增加了 崩塌滑坡点密度这一指标，统计结果见表１． 由表１可看出：汶川Ⅸ度区，地震触发的崩塌 滑坡面积与出现频率之间呈现良好的幂律关系 （Ｒ２ ＝ ０． ９１６）；汶川Ⅹ度区，崩塌滑坡面积与出现 频率之间幂律关系式的相关系数下降到０． ９０６，可 以认为基本服从幂律关系；汶川Ⅺ度区，崩塌滑坡 面积与出现频率之间的关系符合对数正态分布；芦 山Ⅶ度、Ⅷ度、Ⅸ度烈度区，崩塌滑坡规模与发生频 率之间均服从幂律分布，且幂指数ｂ值相近．芦山 地震不同烈度区震害的主要区别是：随地震烈度的 减小，崩塌滑坡点的密度单调减小． 汶川地震、芦山地震是我国迄今为止获得地震 山地灾害资料最全、精度最高的大地震，据此得到 不同烈度区地震触发崩塌滑坡的统计结论，作者更 加关注该结论是汶川地震、芦山地震的独特现象， 还是具有普适性意义的规律．因此，在ＳＯＣ的概念 框架下，利用元胞自动机模拟方法，结合振动台沙 堆模型物理实验，对上述实震资料显示的统计规 律，从物理层面诠释其机理． 图１０ 芦山地震触发的崩塌滑坡分布示意图 Ｆｉｇ． １０ Ｔｈｅ ｄｉｓｔｒｉｂｕｔｉｏｎ ｍａｐ ｏｆ ｌａｎｄｓｌｉｄｅｓ ｉｎｄｕｃｅｄ ｂｙ Ｌｕｓｈａｎ ｅａｒｔｈｑｕａｋｅ ２． ２ 地震触发崩塌滑坡元胞自动机模拟 元胞自动机（ｃｅｌｌｕｌａｒ ａｕｔｏｍａｔａ 或ｃｅｌｌｕｌａｒ ａｕｔｏｍａｔｏｎ，ＣＡ）是空间和时间都离散，物理参量只 取有限数值集的物理系统的理想化模型［１６］ ．该模 型以规则网格形式分布、空间离散的元胞个体为基 本单元，元胞遵循一定的演化规则作同步更新来模 拟真实的物理系统，特别适合复杂系统时空演化过 程的动态模拟研究． 模型方法在ＳＯＣ的研究中具有非常重要的地 位，目前所有对沙堆生长与坍塌机理的模拟算法都 是以元胞自动机为数学基础，对ＳＯＣ 的理解也大 部分来源于元胞自动机沙堆模型的数值模拟．处于 青壮年期的山地系统具有ＳＯＣ的内禀属性，因此， 可以构建地震触发崩塌滑坡沙堆模型，用元胞自动 ９ １ ３
西 南 交 通 大 学 学 报 第５１卷 机来模拟． 根据原型问题的物理特征，与传统沙堆模型相 比，地震触发崩塌滑坡的沙堆模型应具有以下特 点［１７］： （１）传统沙堆模型考察的是一个沙堆在受到 多次扰动下坍塌规模随时间的分布规律．地震发生 时，同时引发多处崩塌滑坡，需要考察多个沙堆在 同时受到一次扰动时坍塌规模的分布规律． 表１ 不同地震烈度区地震触发崩塌滑坡统计 Ｔａｂ． １ Ｓｔａｔｉｓｔｉｃｓ ｏｆ ｌａｎｄｓｌｉｄｅｓ ｉｎ ｚｏｎｅｓ ｗｉｔｈ ｄｉｆｆｅｒｅｎｔ ｓｅｉｓｍｉｃ ｉｎｔｅｎｓｉｔｙ 烈度区崩数塌量滑／个坡崩度塌／（个滑坡／ ｋｍ点２密） 崩塌频滑率坡关规系式模与 检验结果 汶川Ⅺ ３ １５９ — ｆ（Ａ）＝ １ １． ０８ ２槡πＡｅ －（ｌｎＡ － ９． ６０）２ ２ × １． ０８２ χ２值为 １５． １６４，自由度为１０，以χ２ ＜ χ２ ０． ０５（１０）作 为判别标准（查χ２０． ０５（１０）＝ １８． ３０７），认为在显著 水平０． ０５下，服从对数正态分布（９． ６０，１． ０８２） 汶川Ⅹ ２ ８１２ — ｆ（Ａ）＝ ５ × １０７ × Ａ－ １． １３１ 相关系数Ｒ ２ ＝ ０． ９０６，以Ｒ２ ＞ ０． ９为评判标准，认 为服从幂律分布 汶川Ⅸ ３ ３７０ — ｆ（Ａ）＝ １ × １０９× Ａ－ １． ５１３ 相关系数Ｒ ２ ＝ ０． ９１６，以Ｒ２＞ ０． ９为评判标准，认 为服从幂律分布 芦山Ⅸ ４２５ ２． ２５ ｆ（Ａ）＝ ６． ８ × １０４ Ａ－ ０． ９５５ 相关系数Ｒ ２ ＝ ０． ９１６，以Ｒ２ ＞ ０． ９为评判标准，认 为服从幂律分布 芦山Ⅷ ４７７ ０． ４１ ｆ（Ａ）＝ １． ６ × １０５ Ａ－ ０． ９３３ 相关系数Ｒ ２ ＝ ０． ９１７，以Ｒ２ ＞ ０． ９为评判标准，认 为服从幂律分布 芦山Ⅶ ７０６ ０． ３４ ｆ（Ａ）＝ ５． ５ × １０５ Ａ－ １． ０５１ 相关系数Ｒ ２ ＝ ０． ９０１，以Ｒ２ ＞ ０． ９为评判标准，认 为服从幂律分布 （２）在传统沙堆模型中，沙堆是局部受到扰 动，在地震触发崩塌滑坡问题中，坡体是整体受到 扰动，因此，模型需要采用系统整体受扰的方式． （３）传统沙堆模型均是在同一微扰条件下观 察系统的动力学行为．而地震触发崩塌滑坡，首先 不同烈度区，地震对坡体的震动力是不同的，其次 地震力相对坡体已超过微扰量级，因此，设计不同 扰动强度，来模拟这一物理现象． （４）传统沙堆模型，外界对系统的输入是物质 （添加沙粒），发生连锁反应时，扰动的传播是物质 的传播，遵守物质守恒原则．在地震触发崩塌滑坡 的问题中，外界对坡体系统输入地震力，坡体失稳 需要克服其自稳能力，会消耗一部分输入的能量， 扰动的传播是能量的传播，并且存在着能量耗散． 根据以上特点，构造地震触发崩塌滑坡元胞自 动机沙堆模型，基本算法步骤如下： 步骤１ 生成沙堆．地震触发崩塌滑坡相当于 考察大量沙堆，同时受到一次地震作用时坍塌规模 的整体分布规律，因此，一次性生成Ｎ个沙堆，每 个沙堆规模相等，但初始状态不同．对于每个沙堆 而言，考虑一个Ｌ × Ｌ的二维系统，用（ｉ，ｊ）代表元 胞所处的位置（其中，１≤ｉ，ｊ≤Ｌ），每一个元胞有 上、下、左、右４个邻居，Ｆｉ，ｊ为反映元胞（ｉ，ｊ）稳定 性的状态值（相当于元胞的能量）．所有元胞赋予 一个０到阈值Ｆｔｈ之间的初始值，且取随机数． 步骤２ 沙堆演化到临界态．每个沙堆按照下 列规则连续反应，直到所有沙堆演化到临界状态． ① 找到最大状态值Ｆｍａｘ的元胞，把每个元胞 的状态值都增加Δ ＝ Ｆｔｈ － Ｆｍａｘ（相当于对整个系统 的一个扰动），即Ｆｉ，ｊ→Ｆｉ，ｊ ＋ Δ； ② 若新的Ｆｉ，ｊ大于或等于设定的阈值Ｆｔｈ，即 Ｆｉ，ｊ≥Ｆｔｈ，则该元胞倒塌，并向邻居传播扰动，重新 分配Ｆｉ，ｊ给它的４个最近邻： Ｆｉ ± １，ｊ→Ｆｉ ± １，ｊ ＋ αＦｉ，ｊ， Ｆｉ，ｊ ± １→Ｆｉ，ｊ ± １ ＋ αＦｉ，ｊ， Ｆｉ，ｊ→０； （５） ③ 这种扰动传播可能会导致链式反应，如果 由于元胞（ｉ，ｊ）的倒塌导致它的邻居变得不稳定， 重复② 直至所有的元胞状态值小于阈值 （Ｆｉ，ｊ＜ Ｆｔｈ）； 通过步骤２，产生Ｎ 个均处于临界状态，但各 元胞状态值不同的沙堆． 步骤３ 对Ｎ个沙堆同时施加一次扰动．令Ｆ′ 为扰动强度，每个沙堆的元胞状态值都统一增加 Ｆ′，即Ｆｉ，ｊ→Ｆｉ，ｊ ＋ Ｆ′．元胞之间的相互作用仍按照 规则②、③执行． 步骤４ 改变扰动强度Ｆ′，重复步骤３，获得沙 堆在不同扰动强度下坍塌规模与发生频率的关系． ０ ２ ３
第２期 姚令侃，等：山地系统灾变行为自组织临界性研究 需要说明，按照Δ ＝ Ｆｔｈ － Ｆｍａｘ施加扰动时，一 般只会直接触发一两个元胞，但扰动强度加大后， 可能会触发一批元胞，它们各自都可能引发链式反 应，这些链式反应在空间上也会有交叉现象．因此， 在算法上采用并行处理，所有受扰元胞在同一时步 内按照平行更新的方式反应，通过记录倒塌的元胞 数目，作为沙堆坍塌规模的度量． 模型参数取为：α ＝ ０． ２，Ｆｔｈ ＝ １，Ｌ ＝ ５０． 实验生成１００ 万个沙堆（Ｎ ＝ １０６），先连续反 应１０６ 次（均取Ｆ′ ＝ １ － Ｆｍａｘ），以确保沙堆演化到 临界状态．然后取 Ｆ′ ＝ ０． ０００ ０１，０． ０００ ０２，０． ０００ ０４，０． ０００ ０８， １ － Ｆｍａｘ，０． ００１，０． ００５，０． ０１， 开展了８ 组模拟实验．令沙堆坍塌规模为Ｓ，分析 坍塌规模与发生频率的关系，雪崩密度等于发生坍 塌事件的沙堆数除以总沙堆数，实验结果见表２． 表２ 元胞自动机模拟实验结果 Ｔａｂ． ２ Ｒｅｓｕｌｔｓ ｏｆ ｃｅｌｌｕｌａｒ ａｕｔｏｍａｔａ ｓｉｍｕｌａｔｉｏｎ 组 号强扰度动Ｆ′ 发生坍塌事件 的沙堆数／个 总沙堆 数／个 雪崩 密度ρ 坍塌规模 概率关系式 检验结果 １ ０． ０１０ １０６ １０６ １． ００ ｆ（Ｓ）＝ １ ０． ３８ ２槡πＳｅ －（ｌｎＳ － ４． ６８）２ ２ × ０． ３８２ 用χ ２ 检验法在显著水平 ０． ０５下，服 从对数正态分布（４． ６８，０． ３８２） ２ ０． ００５ １０６ １０６ １． ００ ｆ（Ｓ）＝ １ ０． ４５ ２槡πＳｅ －（ｌｎＳ － ３． ８８）２ ２ × ０． ４５２ 用χ ２ 检验法在显著水平 ０． ０５下，服 从对数正态分布（３． ８８，０． ４５２） ３ ０． ００１ １０６ １０６ １． ００ ｆ（Ｓ）＝ １１７． ４１Ｓ－ ２． ９４１ 相关系数Ｒ ２ ＝ ０． ９５１，以Ｒ２＞ ０． ９ 为 评判标准，认为服从幂律分布 ４ １ － Ｆｍａｘ １０６ １０６ １． ００ ｆ（Ｓ）＝ ０． ２６７Ｓ－ ２． ０１３ 相关系数Ｒ ２ ＝ ０． ９６９，以Ｒ２ ＞ ０． ９ 为 评判标准，认为服从幂律分布 ５ ０． ０００ ０８ ４０３ ３０１ １０６ ０． ４０ ｆ（Ｓ）＝ １． ０６８Ｓ－ ２． ４６ 相关系数Ｒ ２ ＝ ０． ９６４，以Ｒ２ ＞ ０． ９ 为 评判标准，认为服从幂律分布 ６ ０． ０００ ０４ ２４２ ５６６ １０６ ０． ２４ ｆ（Ｓ）＝ ０． ５６７Ｓ－ ２． ４４ 相关系数Ｒ ２ ＝ ０． ９６５，以Ｒ２ ＞ ０． ９ 为 评判标准，认为服从幂律分布 ７ ０． ０００ ０２ １４０ ６０７ １０６ ０． １４ ｆ（Ｓ）＝ ０． ３８２Ｓ－ ２． ５０ 相关系数Ｒ ２ ＝ ０． ９６４，以Ｒ２ ＞ ０． ９ 为 评判标准，认为服从幂律分布 ８ ０． ０００ ０１ ８０ ０４９ １０６ ０． ０８ ｆ（Ｓ）＝ ０． ２９９Ｓ－ ２． ５１ 相关系数Ｒ ２ ＝ ０． ９４９，以Ｒ２＞ ０． ９ 为 评判标准，认为服从幂律分布 由表２可看出，以扰动强度１ － Ｆｍａｘ为界可以 划分为两个区间，沙堆模型的动力学特性呈现出不 同的性质：在扰动强度小于１ － Ｆｍａｘ的区间，当Ｆ′从 ０． ０００ ０８ 逐渐减小到０． ０００ ０１ 时，随扰动强度的 降低，坍塌规模与发生频率基本服从同一幂律分布 （表２），但随扰动强度的减小，雪崩密度单调减小 （图１１）；在扰动强度大于１ － Ｆｍａｘ的区间，当Ｆ′从 ０． ００１逐渐增加到０． ０１ 时，随扰动强度的增加，沙 堆模型的动力学特性将经历幂律—幂律弱化—对 数正态分布的渐进式的演变过程（图１２）． 对照两次地震的实震资料统计结论，汶川地震 Ⅸ度区崩塌滑坡规模与出现频率之间的关系可用 幂律描述，Ⅹ度区这一关系基本服从幂律分布， Ⅺ度区更偏向于对数正态．芦山地震ⅦⅨ度烈度 区地震触发崩塌滑坡规模与频率均服从幂律分布， 但随地震烈度的减小，崩塌滑坡点密度单调减小． 通过以上元胞自动机模拟，证明了不同烈度区地震 触发崩塌滑坡的分布规律存在必然的联系性，这些 现象均来源于山地斜坡系统具有ＳＯＣ内禀属性的 物理机制，从而实现了从物理角度对“５·１２”汶川 地震、“４·２０”芦山地震不同烈度区，地震触发崩塌 滑坡分布规律演变机理的诠释，初步建立了不同烈 度区地震触发崩塌滑坡总体特征的完整表征体系． ２． ３ 地震触发崩塌滑坡振动台沙堆模型实验 元胞自动机数值模拟是获得ＳＯＣ系统性质的 主要途径，但相对地震触发崩塌滑坡的原型问题， 无疑是一种高度概化的处理，因此有必要配合沙堆 模型物理实验相互印证． １ ２ ３
西 南 交 通 大 学 学 报 第５１卷 利用大型地震模拟振动台，输入精确可控的地 震波，开展动力扰动下的沙堆模型实验，使之更加 逼近物理原型．实验设备为西南交通大学高速铁路 线路工程教育部重点实验室的电液伺服驱动地震 模拟振动台． 在振动台上放置一长２． ５８ ｍ、宽 １． ５０ ｍ、高１． ９５ ｍ的箱体，选用粒径为０． ６ ～ ５０． ０ ｍｍ经过筛分的干燥天然沙石在箱体里堆成一个 单面坡沙堆，总重量达６． ８ ｔ（图１３）． ρ
图１１ 雪崩密度ρ与扰动强度Ｆ′ Ｆｉｇ． １１ Ａｖａｌａｎｃｈｅ ｄｅｎｓｉｔｙ ρ ｖｓ ｄｉｓｔｕｒｂａｎｃｅ ｉｎｔｅｎｓｉｔｙ Ｆ′ ρ 图１２ 不同扰动强度下坍塌规模的 概率密度曲线 Ｆｉｇ． １２ Ｐｒｏｂａｂｉｌｉｔｙ ｄｅｎｓｉｔｙ ｃｕｒｖｅｓ ｆｏｒ ｄｉｆｆｅｒｅｎｔ ｄｉｓｔｕｒｂａｎｃｅ ｉｎｔｅｎｓｉｔｉｅｓ 图１３ 沙堆模型 Ｆｉｇ． １３ Ｓａｎｄｐｉｌｅ ｍｏｄｅｌ 向振动台输入汶川地震卧龙台站记录的修正 波（图１４），模拟地震对沙堆的扰动．取峰值加速度 为０． ０７５ｇ、０． １００ｇ、０． １２５ｇ、０． １５０ｇ、０． ２００ｇ、 ０． ３５０ｇ、０． ４５０ｇ，开展了７组实验． 图１４ 汶川地震卧龙台记录修正波 Ｆｉｇ． １４ Ｍｏｄｉｆｉｅｄ Ｗｅｎｃｈｕａｎ ａｃｃｅｌｅｒａｔｉｏｎ ｗａｖｅ ｒｅｃｏｒｄｅｄ ｂｙ Ｗｏｌｏｎｇ ｓｔａｔｉｏｎ 每组实验步骤如下： 步骤１ 在箱内堆成单面坡沙堆，坡面靠沙粒 的重力下滑自然形成，使坡脚达到开口端边缘，沙 堆达到临界坡度； 步骤２ 输入指定峰值加速度，称量从坡脚处落 入落沙收集槽的沙粒重量，作为该组的一次实验值； 步骤３ 每次实验后，将落沙收集槽内的沙粒 补充回沙堆上，确保沙堆始终处于临界状态； 步骤４ 重复步骤２、３，每组实验不少于６０次． 以落沙收集槽内沙粒的重量作为规模的度量， 令沙粒重量为ｘ，分析坍塌规模与发生频率之间的 关系，坍塌密度等于发生坍塌事件的实验次数除以 总实验次数，实验结果见表３． 由表３可看出，当输入地震波峰值加速度为 ０． ０７５ｇ、０． １００ｇ、０． １２５ｇ 时，落沙量与发生频率均 服从幂律分布，但随峰值加速度的减小，坍塌密度 单调减小；当峰值加速度增加到０． １５ｇ ～ ０． ２５ｇ时， 落沙量服从对数正态分布；当峰值加速度增加到 ０． ３５ｇ ～ ０． ４５ｇ 时，该阶段样本表现为具有正态分 布的曲线特征． 振动台沙堆模型实验得到了与地震实震现象 规律和元胞自动机模拟类似的结论，更进一步，虽 然目前尚未获得Ⅻ度区的实震资料，但不妨预测应 具有向正态分布发展的趋势． 振动台沙堆模型实验的物理过程更接近地震 触发崩塌滑坡的原型问题，虽然实验次数有限，但 仍然对作者提出的不同烈度区地震触发崩塌滑坡 分布规律的普适性提供了有力支持． ２ ２ ３第２期 姚令侃，等：山地系统灾变行为自组织临界性研究 表３ 振动台沙堆模型实验结果 Ｔａｂ． ３ Ｒｅｓｕｌｔｓ ｏｆ ｔｈｅ ｓｈａｋｉｎｇ ｔａｂｌｅ ｔｅｓｔ ｏｆ ｓａｎｄｐｉｌｅ ｍｏｄｅｌ 组 号加速峰值度／ ｇ 次实数验／次实发验生次坍数塌／的次坍密塌度 坍塌规模 频率 关系式 检验结果 １ ０． ０７５ ９０ ４９ ０． ５４ ｆ（ｘ）＝ ５００． ２ｘ－ ０． ７７４ 相关系数Ｒ ２ ＝ ０． ９００ ５，以Ｒ２ ＞ ０． ９为 评判标准，认为服从幂律分布 ２ ０． １００ ６０ ３５ ０． ５８ ｆ（ｘ）＝ ５７９． ３ｘ－ ０． ７８３ 相关系数Ｒ ２ ＝ ０． ９１７，以Ｒ２＞ ０． ９ 为 评判标准，认为服从幂律分布 ３ ０． １２５ １５０ １１８ ０． ７９ ｆ（ｘ）＝ ３ ８８７． ６ｘ－ １． ０５９ 相关系数Ｒ ２ ＝ ０． ９６２５，以Ｒ２ ＞ ０． ９为 评判标准，认为服从幂律分布 ４ ０． １５０ ６０ ６０ １． ００ ｆ（ｘ）＝ １ ０． ５９ 槡２πｘｅ －（ｌｎｘ － ７． ３６）２ ２ × ０． ５９２ 用χ ２检验法在显著水平 ０． ０５下服从 对数正态分布ＬＮ（７． ３６，０． ５９２） ５ ０． ２５０ ６０ ６０ １． ００ ｆ（ｘ）＝ １ ０． ３２ 槡２πｘｅ －（ｌｎｘ － ８． ３３）２ ２ × ０． ３２２ 用χ ２ 检验法在显著水平 ０． ０５下，服 从对数正态分布ＬＮ（８． ３３，０． ３２２） ６ ０． ３５０ ６０ ６０ １． ００ ｆ（ｘ）＝ １ ５４７１． ５ 槡２πｅ －（ｘ － １６ ０２９）２ ２ × ５ ４７１． ５２ 用χ ２ 检验法在显著水平 ０． ０５下，服 从正态分布Ｎ（１６ ０２９，５ ４７１． ５２） ７ ０． ４５０ ６０ ６０ １． ００ ｆ（ｘ）＝ １ ４ ４８４ 槡２πｅ －（ｘ － １６ ８７９）２ ２ × ４ ４８４２ 用χ ２ 检验法在显著水平 ０． ０５下，服 从正态分布Ｎ（１６ ８７９，４４ ８４２） ３
应
用
案
例
具有ＳＯＣ性质的系统，在临界状态下受到一 系列微小而均匀的扰动，每次扰动下表征反应规模 的物理量可用幂律描述，这种幂律关系将灾害学的 研究与物理学联系在一起．规模与频率的幂次定律 描述的是某种量级灾害发生的数目，这种灾害演变 总体规律在灾势分析、风险评估和危险性区化中的 应用，是ＳＯＣ较为明确的应用领域，对此举一例进 行说明． 除此之外，还将介绍两个属于应用技巧层面的 案例，旨在说明ＳＯＣ存在更广阔的应用前景． ３． １ 岩土体地震扰动深度评估 汶川大地震后，在震区调查时，观察到在无深 大结构面控制的情况下，道路边坡地震失稳一般表 现为浅表的崩塌滑坡，具有“剥皮型”灾害特点．与 日本ＯＹＯ公司合作，采用面波仪对国道２１３ 线都 江堰至映秀段（含水磨支线）震后边坡进行了勘 察，发现自然坡体０ ～ ５ ｍ范围内的表层剪切波速 为２５０ ｍ ／ ｓ，明显低于下部土体的波速，判定为汶 川地震的强扰动范围［１８］ ．但仅以这些表观认识和 有限的测量数据为依据，提出具有普遍性和包容性 的结论是困难的．为此对照１∶ ２ ０００的地形图，对沿 线公路边坡进行详查，内容涵盖沿线公路边坡崩塌 滑坡的规模（崩塌滑坡方量）、崩塌滑坡体深度、工 点经纬度、高程、岩性以及原有工程措施．灾害类型 主要包括砂岩边坡滑坡，岩质边坡崩塌，崩坡积体、 残坡积体、冲洪积体等浅层崩塌滑坡． 调查中发现不同崩塌体的崩塌深度存在很大 区别．在Ⅸ度裂度区调查的６１个崩塌体中，崩塌深 度小于０． ５ ｍ 的２７ 个、０． ５ ～ １ ｍ 的１５ 个、 １ ～ ２ ｍ的１１个、２ ～ ５ ｍ的４个、５ ～ １０ ｍ的３个、 １０ ｍ以上的１个，平均崩塌深度为１． ２ ｍ． 令崩塌滑坡体深度为ｈ，深度大于ｈ的工点数 为Ｎ（ｈ），通过回归分析，拟合曲线如图１５所示．由 图１５可以看出，ｈ与Ｎ（ｈ）在双对数坐标图上近 乎一直线，关系式为 Ｎ（ｈ）＝ ２０． ２７ｈ－ ０． ９９４ ５． （６） 由式（６）可知，相关系数Ｒ２ ＝ ０． ９７３，表明Ⅸ度 裂度区崩塌滑坡深度与出现频率之间存在良好的 幂律关系，说明在Ⅸ度区地震触发崩塌滑坡深度也 存在ＳＯＣ特性．这样若仍以平均崩塌深度作为岩 土体地震扰动深度的表征已不合理． 式（６）为Ⅸ度裂度区崩塌滑坡体深度ｈ与大 于ｈ的工点数Ｎ（ｈ）之间的幂律回归关系，当样本 数量足够大时，在统计意义上可得随机变量ｈ的分 布函数．如据实震调查资料可初步作出ｈ的概率分 布如图１６所示． 由图１６可知，汶川地震Ⅸ度裂度区９３％的崩 塌滑坡体深度均小于５ ｍ，这是基于ＳＯＣ理论开展 ３ ２ ３西 南 交 通 大 学 学 报 第５１卷 的工作，其意义在于建立的是以物理理论为基础的 统计模型，跨越了仅据典型样本得出统计关系的层 面，其规律具有普适性，可以作为汶川地震Ⅸ度区 地震对坡体强扰动效应的定量表征．图１６ 所示的 分布规律可为汶川震后次生山地灾害物源条件估 算以及确定边坡灾害防范深度提供指导． 图１５ 汶川地震Ⅸ度区崩塌滑坡深度累积分布 Ｆｉｇ． １５ Ｌｏｇｌｏｇ ｐｌｏｔ ｏｆ ｌａｎｄｓｌｉｄｅｓ ｎｕｍｂｅｒ ａｎｄ ｄｅｐｔｈ ｉｎ ｚｏｎｅｓ ｗｉｔｈ ｓｅｉｓｍｉｃ ｉｎｔｅｎｓｉｔｙ Ⅸ
图１６ 汶川地震Ⅸ度区崩塌 滑坡体深度概率分布 Ｆｉｇ． １６ Ｐｒｏｂａｂｉｌｉｔｙ ｄｉｓｔｒｉｂｕｔｉｏｎ ｃｕｒｖｅ ｏｆ ｌａｎｄｓｌｉｄｅ ｄｅｐｔｈ ｉｎ ｔｈｅ ｚｏｎｅ ｗｉｔｈ ｓ ｅｉｓｍｉｃ ｉｎｔｅｎｓｉｔｙ Ⅸ ｉｎ Ｗｅｎｃｈｕａｎ ｅａｒｔｈｑｕａｋｅ ３． ２ 泥石流径流量的极值计算 在工程规划与设计中，常常需要对未来工程服 务期内情况进行预测，包括对外荷载最大值的预 计．但人们往往仅掌握有限时段的观察极值，如每 年洪水的最大值．为实现从以前观测到的极值数据 运用外推法求得工程服务期内的极值统计问题，在 数学统计学中已建立了极值统计这一领域，重要成 果之一是来自带有指数型衰减尾部的初始分布的 极值将渐近地收敛于Ⅰ型极限形式［１９］．在以下的 研究中，基于泥石流径流量的频率规模可用幂律 描述关系，通过数学上的映射，将其转换成带有指 数型衰减尾部的分布，从而达到了运用Ⅰ型极限形 式的条件，据此建立了不同设计基准期蒋家沟泥石 流最大径流量极值分布的数学模型． 一场泥石流的总径流量是反映泥石流规模的 主要表征值，也是开展泥石流拦挡工程容量设计、 泥石流堵河可能性预测等工作的关键参数．实际工 程中，关注的是在工程设计基准期内可能出现的泥 石流径流量的极大值，这些极大值可用具有一定概 率分布的随机变量来模拟，但在缺少长序列观测资 料的情况下，确定极值概型是极为困难的． 极值理论是处理一定样本容量极端值分布特 性的理论，其基本概念为：对｛Ｘｔ，ｔ ∈Ｔ｝随机变量 族，将时间域Ｔ等分为ｎ个时段，每个时间段τ ＝ Ｔ ／ ｎ，则得ｎ个独立的随机变量ηｉ（ｉ ＝ １，２，…，ｎ）， 每个随机变量都有相同的分布函数Ｆ（ｘ）．在这个 随机变量族中，极大值为 Ｍ ＝ ｍａｘ（ｘ１，ｘ２，…，ｘｎ）． （７） 在上述两个假定（独立、具有相同的分布函 数）的前提下，极大值的分布函数为 ＦＭ（ｘ）＝ ｐ（Ｍ ≤ ｘ）＝ ｐ（η１ ≤ ｘ，η２ ≤ ｘ，…，ηｎ≤ ｘ）＝ ∏ ｎ ｉ ＝ １ ｐ（ηｉ ≤ ｘ）＝ ［Ｆ（ｘ）］ ｎ， （ ８） 式中：Ｆ（ｘ）为初始分布． 若Ｆ（ｘ）的密度函数为ｆ（ｘ），则有 ｆＭ（ｘ）＝ ｎ［Ｆ（ｘ）］ ｎ － １ ｆ（ｘ）． （９） 由于ＦＭ（ｘ）与初始分布Ｆ（ｘ）有着密切的关 系，这样在随机变量的初始分布以及样本容量ｎ已 知时，可写出极值概率的确切分布．但在通常情况 下，确切分布的求解十分困难，在工程实际应用中， 常用做法就是使用极值的渐近分布． 所谓渐近分 布，就是当数目ｎ趋于无穷大时，概率分布的极限 形式，即渐近分布将收敛于确切分布，这样可用渐 近分布代替实际难以获得的确切分布． 根据Ｇｕｍｂｅｌ的研究［２０］，极值分布的渐近形式 并不取决于初始分布的精确形式，主要决定于初始 分布在极值方向的尾数表现，初始分布的中央部分 对极值分布的渐近形式影响很小，但是，极值参数 则决定于初始分布的形式．工程上常用的渐近分布 有３种类型，其中最大值的Ⅰ型渐近分布多用于荷 载的极值统计方面，如桥梁使用荷载、风速、地震震 级、洪水流量等． 在泥石流研究领域，Ｐ Ａ Ｊｏｈｎｓｏｎ 等研究表 明［２１］，泥石流规模与频率之间存在幂律关系，是自 组织临界状态系统的行为标志．艾南山等认为［２２］ 在泥石流的形成区松散堆积物组元间的非线性作 ４ ２ ３第２期 姚令侃，等：山地系统灾变行为自组织临界性研究 用，使系统自然地朝着临界状态演化，在临界状态 时，诸如降水等外界的细微扰动将被放大而导致规 模不等的泥石流暴发，这种耗散动力学系统的行为 特征，可以用自组织临界状态的概念加以解释．王 裕宜［２３］从泥石流土体的应力应变自组织临界特性 触发，认为粘性泥石流具有自组织临界性．按照 ＳＯＣ系统的定义，泥石流径流量规模Ｗ与频率Ｎ 之间应服从幂律分布，即Ｎ ＝ ａＷ－ ｂ ．假设在所统计 资料中，规模Ｗ的最大值和最小值分别为Ｗｍａｘ和 Ｗｍｉｎ．设Ｘ ＝ ｌｎＷ，则 Ｎ ＝ ａ０ｅ － ｂＸ， （ １０） 式中：ａ０ 为常数；ｂ 为幂指数；Ｘ ≥ Ｘｍｉｎ，Ｘｍｉｎ ＝ ｌｎＷｍｉｎ，一般可取为０． 由频度代替概率的思想，得出Ｘ的分布函数为 Ｆ（ｘ）＝ Ｐ（Ｘ ≤ ｘ）＝ ∫ ｘ Ｘｍｉｎ ａｅ－ｂＸｄＸ ∫∞Ｘｍｉｎ ａｅ－ｂＸ ｄＸ ＝ １ － ｅ－ｂ（ｘ －Ｘｍｉｎ）． （１１） 利用极值分布的最大吸引场的定理［２４］，可以 判定上式极值分布的极限形式为收敛于Ⅰ型的渐 近分布，进而推出Ｘ极大值的极限分布函数为 ＦＭ（ｘ）＝ ｅｘｐ（－ ｅｘｐ（－ α（ｘ － μ）））， （１２） 相应的概率密度函数为 ＦＭ（ｘ）＝ αｅ － α（ｘ － μ）ｅｘｐ（－ ｅｘｐ（－ α（ｘ － μ）））， （１３） 式中：α、μ为一般参数，由如下方法确定：设有ｎ个 单位时间（例如ｎ年）的观测资料，在每个单位时 间内选取一个最大径流量，首先用矩法求得单位时 间最大径流量的平均值珚Ｑ１ 和标准差σＱ １，再根据 式（１４）求得，ｃ２ ＝ １． ２８２ ５５，公式推导见参考文 献［２５］ ． α ＝ ｃ２／ σＱ１ ＝ １． ２８２５５ ／ σＱ１， μ ＝ 珚Ｑｉ－ ０． ５７７ ２２ α ． （１４） 现以蒋家沟泥石流为例，根据蒋家沟泥石流的 实际观测资料（１９６５ ～ １９９４ 年）［２６２８］进行极值分 析．以年为单位，每年选取一个最大径流量，利用矩 法，首先求得年最大径流量的平均值和标准差： 珚 Ｑ１ ＝ ７２８ ４３０ ｍ ３， σＱ１ ＝ ４２１ ６２０ ｍ ３ ．再由式（１４）， 求得参数α和μ，代入式（１２）中，得到蒋家沟泥石 流年最大径流量的概率分布为 Ｆ１（ｘ）＝ ｅｘｐ（－ ｅｘｐ（－ ３． ０４２ × １０ － ６ × （ｘ － ５３８ ６８０）））． （１５） 当设计基准期Ｔ ＝ ５０ ａ时，已知年最大径流量 的分布，统计时段为１ ａ，τ ＝ １． ｎ ＝ Ｔ ／ τ ＝ ５０．根据 极值Ｉ型分布的特征可知，设计基准期最大泥石流 径流量分布的标准差与初始分布相同 σＱ５０ ＝ σＱ１ ＝ ４２１ ６２０ ｍ ３ ． 均值平移距离为 ｌｎ ｎ α ＝ ｌｎ ５０ ／ ３． ０４２ × １０ － ６ ＝ １ ２８６ ００４ ｍ３， 均值为 珚Ｑ５０ ＝ 珚Ｑ１ ＋ １２８ ６００ ４ ＝ ２ ０１４ ４３４ ｍ ３ ． 这样可得设计基准期为５０ ａ 时，蒋家沟泥石 流最大径流量的概率分布为 Ｆ５０（ｘ）＝ ｅｘｐ（－ ｅｘｐ（－ ３． ０４２ × １０ － ６ × （ｘ － １ ８２４ ６８４）））． （１６） 同理可推出设计基准期为１００、３００ ａ时，蒋家 沟泥石流最大径流量的概率分布为 Ｆ１００（ｘ）＝ ｅｘｐ（－ ｅｘｐ（－ ３． ０４２ × １０ － ６ × （ｘ － ２ ０５２ ５４２）））， （１７） Ｆ３００（ｘ）＝ ｅｘｐ（－ ｅｘｐ（－ ３． ０４２ × １０ － ６ × （ｘ － ２ ４１３ ６９０）））． （１８） 概率密度曲线如图１７所示．
图１７ 蒋家沟不同设计基准期的径流量极值概率密度 Ｆｉｇ． １７ Ｐｒｏｂａｂｉｌｉｔｙ ｄｅｎｓｉｔｙ ｃｕｒｖｅ ｏｆ ｅｘｔｒｅｍｅ ｒｕｎｏｆｆ ｖａｌｕｅｓ ｉｎ ｄｉｆｆｅｒｅｎｔ ｄｅｓｉｇｎ ｂａｓｅ ｐｅｒｉｏｄｓ ｏｆ Ｊｉａｎｇｊｉａ ｒａｖｉｎｅ ３． ３ 基于地震活动性参数ｂ值的地应力评估模型 系统动力学服从幂律是人们对ＳＯＣ系统的基 本认识．但幂指数ｂ不仅仅是一个统计分析参数， 还是反映系统动力学特征与活跃度的主要指标，不 同系统ｂ值的差异或者同一系统ｂ值的演化，有着 重要的物理意义．以下基于地震活动性参数ｂ值与 地应力的反相关关系，推求区域地应力的研究，虽 然是一种统计层面的工作，但借此说明，从幂律到ｂ 值的讨论，有可能开拓出ＳＯＣ更广阔的应用空间． 地应力是长大深埋隧道工程设计所需的重要 ５ ２ ３西 南 交 通 大 学 学 报 第５１卷 参数，现场钻探是获得准确初始地应力的可靠手 段，但是在铁路选线的可行性研究阶段，由于线路 方案尚未确定，不宜开展大规模的地应力现场钻探 工作，因此，地应力评估模型的发展就成为具有明 确应用背景的课题． 目前，在地应力场分布特征的研究中多采用统 计实测地应力、建立回归模型的方法来分析主应力 随埋深的分布规律．主要代表性的工作有：中国科 学院武汉岩土力学研究所对中国大陆地区实测地 应力沿埋深分布规律做的相关研究［２９］；中国地震 局地壳应力研究所根据统计的地应力实测数据对 中国陆域八大地块水平主应力随埋深分布规律进 行的回归分析［３０］ ．但是上述研究工作提出的回归 模型，由于公式仅有一个自变参量，即埋深Ｈ，当模 型的统计样本未涵盖评估区域时，就可能造成较大 误差．若能在回归模型中引入能反映当地地应力状 态的参数，则可望通过增加信息量的手段，提高计 算精度． ３． ３． １ ｂ值与地应力的关系 在前述的古登堡里克特定律中，ａ反映平均 地震活动水平，ｂ反映大小地震的比例关系．随后， 在地震预报领域，分析发现，大震前震源及附近区 域经常会出现某些震级档内的地震增多或减小，导 致出现大小地震比例失调、ｂ值减小的异常现象． 由于区域应力积累水平升高是大地震发生的必要 条件，因此，ｂ值也反映了地应力状态，且二者呈反 比关系． 之后有学者［３１３３］根据单轴压缩试验过程中产 生的声发射（ＡＥ）现象对上述假说进行了验证．岩 石的声发射现象描述为：岩石受力变形时，在岩体 内原先存在或新产生的裂缝周围形成应力集中，这 些局部应力集中区不均匀发展，并发生突然的破 裂，从而向四周辐射弹性波．声发射活动与地震活 动的机制最为接近，在统计参数上与地震活动性的 可对比性也最强，具体表现为：实验记录到的声发 射次数与声压幅值关系服从幂律分布，其中，幂律 关系式中的幂指数ｂ值不仅为一常数，而且随着应 力的增加明显下降，即岩石处于高应力状态时， ｂ值低．例如从中国国家地震局地球物理研究所公 布的实验资料（图１８）中可看出，岩石破裂实验声 发射ｂ值随应力的增加而降低，在达到破裂应力之 前ｂ值下降较快［３４］． 基于ｂ值与地应力相关的原理，由于ｂ值是可 以通过历史地震资料统计获得的参数，因此，在水 平主应力随埋深分布规律统计回归方程的基础上， 可以增加ｂ值作为反映各地地应力状态的信息参 量，对现有地应力评估模型进行改进．
σσ 图１８ ５种岩石声发射ｂ值随应力变化图［３４］ Ｆｉｇ． １８ Ｃｈａｎｇｅ ｗｉｔｈ ｓｔｒｅｓｓ ｏｆｔｈｅ ｂｖａｌｕｅｓ ｏｆ ａｃｏｕｓｔｉｃ ｅｍｉｓｓｉｏｎ ｆｒｏｍ ｆｉｖｅ ｔｙｐｅｓ ｏｆ ｒｏｃｋ［３４］ ３． ３． ２ 基于ｂ值的川藏交通廊道地应力评估模型 川藏交通廊道穿越印度洋板块俯冲欧亚板块 形成的青藏高原，板块边界的作用力是构造变形的 主要动力源，同时造成板块内部地震活跃、大地变 形、高地应力等现象［３５］ ．如图１９所示，从青藏高原 东缘进藏，需穿越三大缝合带（金沙江缝合带、班 公怒江蛇绿岩缝合带、雅鲁藏布江蛇绿岩缝合带） 和七大断裂带（龙门山断裂带、鲜水河断裂带、甘 孜理塘断裂带、金沙江红河断裂带、甲桑卡赤布 张错断裂带、纳木错仲沙断裂带、嘉黎至然乌断裂 带），构造活动较强烈．此外，进藏线路靠近龙门山 地震带，穿越甘孜炉霍地震带、雅鲁藏布江地震带， 地震烈度均在Ⅷ度以上，地震活动强烈． 以北纬２８° ～ ３２． １°、东经９０． １° ～ １０４． １°范围 内，１９７０年１月～ ２０１３年１２月共４５ ａ的地震资料 作为基础数据，震中分布见图２０．将地震记录以 ０． ２ 震级档为间隔进行分级整理．利用Ａｒｃｇｉｓ软件 将研究区平面图以０． １° × ０． １°的间距进行网格 化，挑选出以每个网格节点为圆心、半径为ｒ的圆 形统计单元内的地震资料［３６］，对于多数节点，统计 单元的半径ｒ ＝ ３０ ｋｍ，但由于地震活动程度存在 较大差异，对于地震分布较稀疏的局部区域，统计 单元ｒ值可适当扩大，以满足统计所需的地震样本 数．对于每一个统计单元，采用非线性拟合方法计 算出拟合公式以及拟合优度Ｒ２，确定能满足整个 研究时段的最小完整性震级Ｍｃ［ ３７］，不同统计单元 的Ｍｃ 是有差异的，计算中选择拟合优度Ｒ ２ 取最 大值时对应的震级为Ｍｃ．利用最佳拟合优度法，得 到最小完整性震级取Ｍｃ 时的地震活动性参数ｂ ６ ２ ３
第２期 姚令侃，等：山地系统灾变行为自组织临界性研究 值，作为该单元中心点（即网格节点）的计算ｂ值． 利用Ａｒｃｇｉｓ软件将得到的整个研究区域扫描计算 网格点的ｂ值对离散点进行插值计算，经数值范围 和区域颜色的调整后得到ｂ值的空间分布等值线 图（图２０）．由图２０ 可知，ｂ值空间分布沿构造断 裂带存在着明显的空间差异，反映出不同断裂带以 及同一断裂带不同断裂段现今应力积累水平的差 异．川藏交通廊道中，雅安、八美、道孚、甘孜、巴塘 至二九六工班、邦达至八宿、通麦、墨竹工卡、曲松 等地为异常低ｂ值区，地应力积累水平较高． 图１９ 川藏交通廊道线路、断裂带、缝合带、中小地震震中分布图 Ｆｉｇ． １９ Ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ ｌｉｎｅｓ，ｆａｕｌｔｓ，ｓｕｔｕｒｅｓ，ａｎｄ ｅａｒｔｈｑｕａｋｅ ｅｐｉｃｅｎｔｅｒｓ ｏｎ ｔｈｅ ｔｒａｎｓｐｏｒｔａｔｉｏｎ ｃｏｒｒｉｄｏｒ ｆｒｏｍ Ｓｉｃｈｕａｎ ｔｏ Ｔｉｂｅｔ 图２０ 川藏交通廊道ｂ值空间分布 Ｆｉｇ． ２０ Ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ ｂｖａｌｕｅｓ ｏｆ ｔｈｅ ｔｒａｎｓｐｏｒｔａｔｉｏｎ ｃｏｒｒｉｄｏｒ ｆｒｏｍ Ｓｉｃｈｕａｎ ｔｏ Ｔｉｂｅｔ 本文中收集了１０ 个测区，１０２ 组地应力原始 实测数据．所收集的原始地应力数据测量方法多为 钻孔应力解除法和水压致裂法，少数为凯瑟尔效应 法．地应力测试工程主要包括：南水北调西线一期 工程、二郎山隧道、嘎隆拉隧道、大岗山水电站、两 河口水电站等．各测点所在区域的ｂ值见表４． 利用上述资料进行统计回归分析，建立了实测 点最大主应力σ１、埋深Ｈ、地震ｂ值三者的关系 式，即川藏交通廊道地应力评估模型为 σ１ ＝ ２． ２７０ ７６ － １． ４１７ｂ ＋ ０． ５００ ４Ｈ． （１９） 表４ 实测点与ｂ值统计表 Ｔａｂ． ４ Ｓｔａｔｉｓｔｉｃｓ ｏｆ ｍｅａｓｕｒｉｎｇ ｐｏｉｎｔｓ ａｎｄ ｂｖａｌｕｅｓ 图中编号 １ ２ ３ ４ ５ ６ ７ ８、９、１０ １１ １２ 测点 冈南底ＣＫ３斯地钻块孔 嘎隧隆道拉ＺＫ０９泥曲钻孔甘孜石绒门岔坎寺、水两河电口站硗碛宝兴 康定 二隧郎道山大水岗电山站 ｂ值 ２． ５ ２． ５ １． ５ ０． ６１１ １ １． ０３５ ０． ６３ １． ０４ １． ０３ ０． ８９ ７ ２ ３
