2020, Vol. 63
2. 中国矿业大学 安全工程学院, 江苏徐州 221116;
3. 中国矿业大学 煤矿瓦斯与火灾防治教育部重点实验室, 江苏徐州 221116
2. School of Safety Engineering, China University of Mining and Technology, Xuzhou Jiangsu 221116, China;
3. Key Laboratory of Gas and Fire Control for Coal Mines(China University of Mining and Technology), Ministry of Education, Xuzhou Jiangsu 221116, China
岩体内部孕含着微宏观多尺度多层次的裂隙，这些岩体结构在载荷作用下裂纹扩展行为及力学响应关系着地下工程的稳定性及其失稳致灾过程(Lennartz-Sassinek et al., 2014; Moradian et al., 2016).然而裂隙岩体裂纹扩展是一种复杂的断裂力学行为，由于裂纹尖端应力场的奇异性，对单裂纹在理想介质中的起裂、扩展问题也无法用线弹性断裂力学完美解释(Nanjo, 2017).因此，非线性断裂力学中的塑性区概念、Dugdale-Barenblatt模型、等效弹性裂纹等模型被提出并广泛使用(Lin and Labuz, 2013).载荷作用下裂隙岩体裂纹起裂、扩展、成核等行为不仅仅取决于宏观应力场及损伤场的分布性质，还很大程度上受到微观尺度破裂活动的影响(Graham et al., 2010)，这就要求我们必须准确地掌握这微破裂的空间位置、张剪属性、几何尺度、空间取向及运动方向等微观属性.
在实验室尺度下，预制裂纹岩石受载微破裂活动属性的刻画和反演可以借助声发射(acoustic emission, AE)监测来实现(王笑然等，2018).岩石微破裂处产生位错扰动，该扰动以震动波形式向四周传播，被介质表面的传感器检测到，记录的扰动响应即为声发射波形信号(Ishida et al., 2017; Wang et al., 2017).声发射信号的产生机理、传播方式、波形特征及震级频度幂次法则都与大尺度地震波类似，因此微破裂活动也常被认为“微缩版”的地震(刘培洵等，2014).由于声发射与地震波的相似性，定量地震分析法常被借鉴到声发射信号处理中，形成了一整套定量声发射分析理论(Rice, 1980; Othsu, 1991, 1995; Shah and Labuz, 1995；Carvalho and Labuz, 2002；Kao et al., 2011;刘培洵等，2014).声发射定量反演属于弹性动力学范畴，即通过提取各传感器处波形数据并结合定位结果能获取破裂源类型及产生机制等一系列震源参数.然而，声发射原始波形数据不仅包含微破裂震源参数信息，同时也包含了波在介质中的传播及接收系统的响应信息.为了获得微破裂的震源参量，必须将震源模型信息从波传播和传感器响应信息中定量分离，即通常需要将声发射定量反演研究细分为震源模型、波在介质中传播、传感器响应三个具体问题(Carvalho and Labuz, 2002).
震源机制定量反演首先要建立破裂源产生模型.在岩石等准脆性材料中，声发射源通常是由于材料体内出现了位移不连续(displacement discontinuity, DD)或应变不连续(strain discontinuity, SD).位移不连续模型是仅考虑裂纹张开/闭合和剪切滑移的微裂纹震源模型(Kao et al., 2011)，有些学者也称其为剪切-张拉裂纹(shear-tensile crack, STC)模型(Petružálek et al, 2018)；而应变不连续模型用来描述微体积等向变化产生的外爆或者内塌震源(Shah and Labuz, 1995).在定量地震学中这两种模型被等同对待，即位移和应变不连续都可等效为一系列力偶的叠加(Aki and Richards, 1980).对于声发射在介质中传播，通常将其简化为无限大线弹性各向同性均质材料中埋入式扰动脉冲在空间的位移场响应问题，这样就可以采用现有的格林函数进行求解：Stokes解给出了单极体力造成的动态格林函数，在此基础上Aki和Richards(1980)计算得到了由矩张量造成位移场的双力偶解，Othsu(1991)编写了专门针对声发射矩张量反演的简化格林函数分析SIGMA程序并取得了广泛应用，这些研究为利用矩张量反演震源机制奠定了理论基础.通过对比传感器位置处的实测P波幅值数据和矩张量求解的理论位移数据，可以采用最小二乘法求解震源矩张量分量，在这种情况下，传感器响应及标定成为震源反演的重要环节.
基于上述三方面综合考虑，大量声发射震源机制定量反演理论问世，并在岩体失稳监测、水力压裂机理研究等方面取得了很多成果(Othsu, 1991; Shah and Labuz, 1995；Carvalho and Labuz, 2002；Chang and Lee, 2004; Liu et al., 2015, 2018; Wong and Xiong, 2018; Hampton et al., 2018).但目前常用的震源机制反演程序也存在着一些问题与不足：(1)震源模型的建立与分解前后矛盾或物理意义不明确.如Othsu(1991, 1995)虽从位移不连续模型出发，但其将矩张量特征值分解为各向同性(ISO)、纯双力偶(DC)和补偿线性矢量偶极(CLVD)三部分，然后根据各部分占比来判断微破裂的属性.DC部分的物理意义可解释为纯剪切破裂，但ISO和CLVD部分在震源模型中的物理意义不明确，当且仅当ISO和CLVD按特定占比组合才能将其物理意义解释为纯张拉破裂.为此，Shah和Labuz(1995)将矩张量特征值分解为纯剪切(DC)部分和纯张拉(ISO和CLVD按特定比例组合)部分，这两部分组合就是剪切-张拉微裂纹震源模型(位移不连续模型)；而余下的ISO部分可解释为微体积等向变化产生的外爆或者内塌震源(应变不连续模型).即他们认为震源是微裂纹事件和微体积改变的线性叠加，很显然，这种处理并不太令人满意，事实上震源通常只能是位移或应变不连续两者中的其一，而不是两者的叠加组合.特别是对于低孔隙硬岩，裂纹扩展产生的声发射事件通常是位移不连续造成的，这就需要我们建立考虑裂纹张开/闭合和剪切滑移的位移不连续震源模型.(2)缺少对各位置处传感器耦合质量的定量评价和标定.传感器响应是震源反演的重要环节，而传感器采集到的数据很大程度上取决于传感器与试样界面之间的耦合质量，为了确保各传感器与试样的耦合质量一致，通常采用相对标定实验(刘培洵等，2014)，通过相对位移来获取相对矩张量而无法得到其绝对大小，故无法与全局监测的其他参量(如应力、应变等)对比验证.此外，传感器采集数据通常为电压值非位移量，即传感器响应是转换方程与该位置处扰动的卷积，通常采用去卷积法复原传感器处的位移数据，但该过程计算复杂耗时，且大部分为无用功，因为我们感兴趣的仅是P波起跳幅值.
基于此，本文以含预制裂纹砂岩试样为对象，测试其在单轴压缩载荷下宏观裂纹扩展及伴随的声发射信号，采用标量敏感参数法对声发射传感器耦合质量进行绝对标定，建立仅涉及微裂纹张开闭合和剪切滑移的位移不连续震源模型，通过点源远场P波矩张量反演获得微裂纹类型、体积、空间取向及运动方向等微观参量，分析裂纹扩展过程中这些微观参量的时变响应、统计特征及微裂纹活动的断裂力学行为.研究成果可为进一步认识裂隙岩体裂纹扩展微观机制提供参考.
1 实验方案 1.1 试样及性质本文实验试样为红砂岩，试样致密且均质性好.通过电动切割机沿同一方向切割加工成长方体试样，用磨石机反复磨平以确保试样端部不平整度小于0.2 mm；然后用高压水射流预制两条平行贯穿裂纹，制备好的砂岩试样尺寸见表 1和图 1.采用波速测试系统(ARB-1410卡、WaveGen控制软件、PCI-Ⅱ波形显示软件)通过透射脉冲法测试并计算试样在x, y, z三个方向的P波和S波波速(Wang et al., 2019)，并通过式1计算试样各向异性系数ΔA(Cerrillo et al., 2014)，结果均在3.5%以内(表 2)，故可以认为试样为各向同性均质材料.根据砂岩试样P波波速cP、S波波速cS及密度ρ，计算其弹性模量EY、拉梅常数λ和μ和泊松比ν，结果如表 2所示.
|
图 1 试样及实验系统图 (a)试样示意图；(b)预制裂纹几何参数；(c)测试系统. Fig. 1 Diagram of the specimen and experimental system (a) Sandstone specimen; (b) Geometric of the pre-existing cracks; (c) Test system. |
|
表 1 砂岩试样尺寸及声发射传感器位置坐标 Table 1 The dimensions of sandstone and coordinates of AE sensors |
|
|
表 2 砂岩试样的弹性常数 Table 2 Elastic constants of the chosen sandstones |
|
(1) |
其中，UPV(Ultrasonics Pulse Velocity)为超声脉冲波速，把P波波速cP或S波波速cS代入可计算由P波或S波得到的各向异性系数ΔA(见表 2)；下标i, j, k代表x, y, z三个方向且i≠j≠k.
1.2 实验系统及方案本文测试系统主要包括加载系统和声发射信号采集系统.加载系统为MTS伺服压力机(最大载荷600kN)，实验中采用常位移速率(0.15 mm/min)加载控制；声发射系统包括24通道美国物理声学公司Micro-Ⅱ型采集主机、前置放大器和NANO-30声发射探头，能实时采集试样受载中声发射原始波形数据，并通过编写Matlab程序对声发射事件进行P波到时及起跳幅值自动拾取、3D定位及震源机制反演.在试样前后表面不同位置各布置4个声发射传感器(空间坐标见表 1)，来确保声发射事件尽可能地落入传感器阵列中.实验开始前采用断铅试验测试传感器耦合质量及监测系统精度，设置前置放大器增益为40 dB，触发阈值为40 dB，信号采样率为1 MHz，采样长度为1024个点(其中前1/4为预触发).加载过程中，实时同步采集含预制裂纹砂岩试样的力学(载荷、位移等)和声发射波形数据，并同时采用Canon EOS 6D相机(采样率24 fps)对试样表面裂纹扩展过程进行录像.
2 声发射到时拾取及定位算法 2.1 P波到时及起跳幅值自动拾取声发射震源机制定量反演的基础是可靠的声发射定位数据，而高精度定位结果获取的前提是精确的到时拾取.声发射到时拾取就是确定P波震相的初到时刻(Allen, 1982)，而到时之前的信号被认为是背景噪声，在此基础上各种各样的P波到时自动拾取算法(Maji and Shah, 1988; Kurz et al., 2005; Carpinteri et al., 2012;白添羊等，2016)被提出.其中常用的到时拾取算法有：①噪声阈值法，即计算背景噪声信号的平均值、方差等统计量来确定门槛阈值，基于声发射信号第一次越过该阈值的时刻向前线性外推来确定P波到时；②长短时窗法(SLA/LTA)，即计算信号特征包络函数的SLA(Short Term Averaging)和LTA(Long Term Averaging)值，然后由SLA与LTA比值来确定P波到时；③赤迟信息准则法(Akaike Information Criterion，AIC)，即通过自回归模型计算选用时窗内信号AIC函数值，AIC最小值处为背景噪声信号与P波信号的分界点，该点即为到时点.为了尽可能精确地确定P波到时，本文以AIC法为基础，到时拾取过程主要包含三个步骤：
(1) 背景噪声信号统计量及阈值确定：选取一段背景噪声信号Sn(如图 2a中前200个点)，计算Sn的平均值m以及标准差σ，则噪声信号阈值Tn为
|
图 2 声发射到时确定及其不同方法对比 (a)声发射P波到时确定步骤；(b)真实背景噪声；(c)背景噪声放大3倍. Fig. 2 First arrival determination based on the proposed AIC approach and comparison with other methods (a) Procedure for determining P-wave arrival time; (b) Situation of actual environmental noise level; (c) Situation of three times environmental noise level. |
|
(2) |
其中：Ln为选取的背景噪声Sn的长度.
(2) 基于AIC法确定P波到时：记声发射波形首次越过噪声阈值Tn的点为Kf, 以Kf为目标点向前选取Kb个点，向后选取Ka个点，由式3计算时间窗[Kf-Kb, Kf+Ka]内声发射片段R(图 2a)的AIC(w)函数值.Ka和Kb值的选取须确保片段R中既包含部分背景噪声信号又包含P波信号，也就是说P波到时点必须包含在R中.记AIC(w)函数值在wmin处取得最小值，则潜在P波到时点为wmin+Kf-Kb-1.
|
(3) |
其中，LR为选取信号片段R的长度，LR=Ka+Kb+1; w为选用片段R的索引号，w取值范围为[1, LR]; var为样本方差函数(var=σ2).
(3) 到时点验证：为了确保步骤(2)中的潜在到时点为真实P波到时点，需同时验证①从该点开始在之后30个点内，声发射幅值至少跨越4次背景噪声信号阈值Tn；②从该点开始在之后50个点内，声发射幅值至少改变符号5次.当同时满足上述两个条件，则根据该到时点、采样率及触发时刻，确定并记录P波到时；否则，该到时数据被舍弃.
为了说明该方法优越性，在相同P波信号数据、不同背景噪声下，将该法拾取的P波到时与噪声阈值法、长短时窗法及手动拾取进行对比，结果如图 2b和2c.可以发现当噪声信号幅值较低时，这四种到时拾取法获得的P波到时都接近或者完全一致.而当背景噪声信号幅值增大三倍后，噪声阈值法和长短时窗法与手动P波到时拾取结果相差较远，而本文的AIC法仍能得到较准确的P波到时.
P波到时确定后其起跳幅值便可随之确定，即满足式(4)的第一个信号幅值Ri(i=wmin, …, LR)为P波幅值VP.根据以上步骤，编写Matlab程序对声发射P波到时及起跳幅值进行自动拾取(结果见图 3)并存储，其中P波到时数据为声发射定位的输入文件，而P波起跳幅值数据为震源机制反演输入文件.
|
图 3 典型声发射波形及其P波到时与起跳幅值确定 Fig. 3 The time histories of one AE event and P-wave onset times and amplitudes |
|
(4) |
声发射源尺度远小于震源到各传感器间的距离，故假设其为点源，则声发射源定位控制方程为
|
(5) |
其中，x (x, y, z)和xi(xi, yi, zi)分别为声发射源及第i个传感器的空间坐标；cP为P波波速；ti为第i个传感器处P波到时；t为声发射源的发震时刻.
当P波波速已知时，式(5)中含有(x, y, z, t)四个未知数，考虑到波速模型、P波到时拾取不准等随机因素，通常采用尽可能多的传感器建立超定方程组，采用最小二乘法来消除随机因素(Ge, 2003).该方法在随机误差(近似)服从高斯分布时效果较好，而当采集数据中存在大误差离群点时，定位精度将大打折扣.因此，最小绝对偏差(又称L1范数统计)定位算法受到重视(刘培洵等，2009；李楠等，2014)，该算法把各个传感器处采集到时与计算到时之差绝对值和作为目标函数(式(6))，使目标函数f(x)取得最小值的点为震源坐标解：
|
(6) |
在最小绝对偏差法中，震源发震时刻t的最佳估计为各个传感器处采集到时与计算到时之差的中位数(Ge, 2003).本文采用单纯形数值优化算法搜索目标函数的最优值，即计算“单纯形”各个顶点处的f(x)值，确定最优点B(f(x)值最小)及最差点W(f(x)值最大)，计算除W点之外其他顶点构成图形的形心O；通过基本变形运算得到W点关于O点的映射点R、外扩点E、内缩点C，比较W点与R、E、C点处f(x)值的大小，如果W点为R、E、C点中的最优点，则将原单纯形朝B点等比收缩(S)形成新的单纯形，否则将W点替换成R、E、C中的最优点，组成新的单纯形.重复上述步骤，直到单纯形大小满足迭代终止条件，此时单纯形最优顶点为目标函数的最优解.单纯形在解空间中总是朝着函数值最小的方向移动，故不会出现迭代不收敛问题；此外，单纯形搜索法虽迭代步数较多，但其避免了偏导数和逆矩阵的求解，既降低了计算量又防止了Hessian矩阵不可逆等问题，使得计算过程更加快速稳健.
根据声发射P波到时数据，采用单纯形定位稳健算法对声发射源进行定位.在FS-0.15-2试样实验中共记录到6455个声发射事件，其中2198个事件被定位.为了定量评价定位事件的精度，本文定义并计算各个事件定位误差(式(7))，其分布如图 5e所示.其中90%以上定位事件的误差在5 mm以内，1149个(52.27%)事件定位精度在2.5 mm以内(其空间分布见图 5(a—d)).这些定位事件密集区与裂纹扩展路径、试样破坏宏观主裂纹(图 5f)相吻合，这表明了本文到时拾取及定位算法的可靠性.
|
图 4 单纯形搜索法及其四种基本变化(Ge, 2003) Fig. 4 An illustration of the simplex method and the four basic deformation mechanisms (Ge, 2003) |
|
图 5 FS-0.15-2试样声发射定位及破坏图 (a) 3D声发射定位结果; (b) X-Y平面投影; (c) Z-Y平面投影; (d) X-Z面投影; (e)定位误差分布; (f)试样最终破坏模式. Fig. 5 AE source locations and failure mode of FS-0.15-2 specimen (a) 3D AE source location results; (b) X-Y plane projection; (c) Z-Y plane projection; (d) X-Z plane projection; (e) AE source location error distribution; and (f) Final failure mode. |
|
(7) |
为了表征低孔隙硬岩裂纹扩展震源机制，本节建立仅涉及裂纹张开/闭合和剪切滑移的位移不连续震源模型.设产生声发射的微裂纹事件表面为Γ，其法向为n，由微裂纹位错产生的位移不连续向量(Burgers向量)为b，则该微裂纹运动学特征可以n和b定量描述.如果把震源当作点源处理，则由该微裂纹事件产生的地震中常用的矩张量为(Othsu, 1995)
|
(8) |
其中，ΔA为位移不连续面积；Cijpq为材料弹性刚度张量.对于各向同性均质材料，Cijpq可表示为
|
(9) |
其中，λ和μ分别为材料的拉梅常数；δij为Kronecker delta函数.
把式(9)代入式(8)中，则各向同性均质材料内由位移不连续产生的矩张量为
|
(10) |
其中，k为哑标，即bknk=b1n1+ b2n2+ b3n3.
由于矩张量Mpq为二阶对称张量(式(10))，其独立分量可减少至6个.由特征值和特征向量定义：
|
(11) |
其中，M为矩张量的特征值(主值)；xi为特征向量(主轴方向)，i为自由下标.
式(11)有解的充要条件为(Mpq-Mδpq)的行列式值为0.对其行列式展开可得到关于M的3次方程，方程的3个根即为矩张量的3个特征值，回代到式(11)中即可求得矩张量的三个特征向量，结果为
|
(12) |
其中，M(1)、M(2)和M(3)分别为最大、中间和最小特征值(M(1) ≥M(2) ≥M(3))；xi(1)、xi(2)和xi(3)分别为对应的特征向量；|b|为位移不连续向量的模，即|b|= 
b和n夹角(裂纹模式角)α决定了微裂纹的类型.当α=0°或180°时，b//n，微裂纹为Ⅰ型张开或闭合裂纹；当α=90°时，b⊥n，微裂纹为纯剪切(Ⅱ型滑开或Ⅲ型撕开)型裂纹；当0° < α < 90°时，微裂纹为张开和剪切混合裂纹；而当90° < α < 180°时，微裂纹为闭合和剪切混合裂纹.由式(12)可得裂纹模式角α和微裂纹体积Vol与矩张量特征值之间的关系，即
|
(13) |
|
(14) |
由xi(2)=εijkbjnk可知特征向量x(2)与b和n所在平面垂直，且特征向量x(1)平分b和n的夹角(见图 6).这里需要注意的是，b和n空间取向并无法唯一确定(图 6(a, b)两种情况)，即b和n互换可以产生相同的矩张量(式(10))，因此需要提供额外的信息来唯一确定b和n，这部分内容在第4节中进行讨论.
|
图 6 矩张量主轴坐标系下微裂纹法向和位移不连续向量取向 (a)情况1；(b)情况2. Fig. 6 Normal and displacement discontinuity vectors under eigenvector coordinate of the moment tensor (a) Case 1; and (b) Case 2. |
需要注意的是在仅考虑裂纹张开/闭合和剪切滑移的位移不连续震源模型，由式(12)可知矩张量的3个特征值与弹性常数之间必须满足某种约束条件(Carvalho and Labuz, 2002; Kao et al., 2011)，即：
|
(15) |
其中，ν为试样的泊松比.
3.2 点源远场P波矩张量反演为避免复杂表面反射波及S波等后续震相的影响(刘培洵等，2014)，本文采用P波对矩张量进行反演.根据定量地震学理论，在无限大各项同性均匀介质中由点源P波矩张量Mpq(x 0, t0)在任意位置x处、在t时刻产生的远场径向位移为(Aki and Richards, 1980)
|
(16) |
其中，Grp, q为格林函数的空间导数；*为卷积运算符；x0为声发射源位置向量；t0为声发射源发震时刻，t0=t-r/cP; γ为声发射源到传感器向量的方向余弦；ρ为介质密度；cP为P波波速；r为声发射源与传感器的空间距离；
|
(17) |
在式(16)，仅P波远场项位移被考虑，该近似是合理的.对于本文试样尺寸(65 mm×30 mm×150 mm)，声发射源与传感器之间的距离通常大于几倍声发射波长(波速/主频=5000[m/s]/300[kHz]≈16.7 mm).传感器处采集的P波波形取决于矩张量随时间的变化率
|
图 7 传感器耦合质量绝对标定测试原理 (a)基于断铅测试的法向位移计算示意图; (b)传感器采集的波形及理论法向位移对比. Fig. 7 Principle of absolute calibration test for sensor coupling quality (a) Normal displacement calculation based on lead-break test; (b) Comparison between calculated normal displacement and measured AE signal. |
|
(18) |
其中，Tu为微破裂位错时间.Carvalho(1999)通过Lamb弹性理论解和传感器标定测试获得了Tu与介质弹性常数的函数关系：
|
(19) |
由于P波传播方向与质点粒子的运动方向一致，故认为P波传感器仅对法向垂直位移分量产生响应.而P波传播到试样与传感器交界面时，由于入射P波与传感器法向存在一定的夹角，会产生反射P波及由P波模式转变成的SV波(Scruby et al., 1985)，则径向入射P波幅值ur与传感器采集波形法向幅值un的关系为
|
(20) |
其中，Rp为反射系数；β为入射P波与传感器法向夹角；k为P波与S波波速比.
3.3 传感器耦合质量标定声发射传感器采集到的P波幅值数据很大程度上取决于各个位置处传感器与试样间的耦合质量.为了消除传感器耦合质量对实验结果的影响并获取矩张量的绝对大小，且避免复杂繁琐的去卷积过程，本文采用标量敏感参数标定实验法(Carvalho, 1999)，即实验开始前在试样表面不同位置进行断铅测试，然后通过式(21)计算(原理见图 7a)任意泊松比下Lamb问题的数值解来获得各个传感器位置处的理论法向位移un(τ)(Mooney, 1974).un(τ)的理论结果为图 7b中的波形虚线，其第一个峰值处的位移up为P波理论法向位移.
|
(21) |
其中，P为铅芯断裂所需的力; τ为无量纲时间数，τ=cSt/r；r为断铅点与传感器间的距离；G(τ)为理论法向位移的时间依赖函数(Mooney, 1974).
采用断铅试验中传感器记录的P波起跳幅值VP的绝对值及理论计算幅值up，定义幅值标量敏感参数SA为
|
(22) |
实验中实际采集的声发射P波起跳幅值Vp乘以SA便可获得传感器处采集到的绝对法向位移值.通过SA既定量化了传感器与试样间的耦合质量，也代替了复杂的去卷积过程，使得声发射震源定量反演过程大大简化.
3.4 模型求解及流程把式(17)—(20)代入式(16)中，可得任意位置处在t=Tu+r/cP时刻下由矩张量计算的理论法向位移.求解矩张量的6个未知数分量，使得各位置处传感器采集的P波起跳幅值与该位置处矩张量计算位移之差平方总和最小：
|
(23) |
其中，(SAVp)i和(urRp)i分别为第i个传感器位置处采集的法向位移和矩张量计算的法向位移.
在3.1节中建立的仅考虑裂纹张开/闭合和剪切滑移的位移不连续震源模型中，必须满足等式(15)的约束条件，即位移不连续模型需求解在等式(15)约束条件下等式(23)的最小值问题.考虑到式(15)的约束条件是在矩张量的主轴坐标下，而矩张量分量求解是在任意坐标系下，因此需要把式(15)的约束条件改写为矩张量不变量的形式(与应力张量不变量类似)：
|
(24) |
其中，I1、I2、I3分别为矩张量第一、第二和第三不变量，其值与坐标系的选取无关.I1、I2、I3与矩张量的各个分量关系为
|
(25) |
本文采用拉格朗日乘子法求解该约束条件极值问题，即通过引入新的未知数乘子γ后转化为求解新目标函数EN=E+γEC的最小值问题.虽然原始目标函数(式(23))为典型的凸函数优化问题，但由于高阶非凸约束条件的存在(式(24))，使得新的目标函数求解变得困难.这里采用Levenberg-Marquardt迭代算法求解该约束条件优化问题，其核心迭代更新公式为
|
(26) |
其中，xk和xk+1分别为第k步和k+1步的未知数向量(x=[x1 x2 x3 x4 x5 x6 x7]T, x1~x6为矩张量的6个未知数分量，x7为拉格朗日乘子γ)；g为新目标函数EN的梯度向量，gi=∂EN/∂xi(i=1 : 7)；H为新目标函数EN的海森矩阵，Hij=∂2EN/∂xi∂xj (i, j=1 : 7);ξ为阻尼调整因子，在每步的迭代过程中根据EN值的增减来调整ξ的大小：当ξ很大时，迭代公式趋近于x k+1= x k- g /ξ, 为梯度下降法；当ξ很小时，迭代公式趋近于xk+1= xk- H-1 g，为牛顿法.
Levenberg-Marquardt迭代算法通过引入阻尼因子ξ既提高了迭代的收敛速度又有效地解决了Hessian矩阵不可逆问题.当梯度向量g的L2范数小于设定的迭代收敛精度(1×10-10)时，迭代终止，此时的x1~x6即为矩张量的6个分量，进而根据式(12)—(14)求解微裂纹类型及产生机制等震源参量.
综合第2节和第3节内容，本文编写了Matlab程序对砂岩试样受载过程中的震源机制进行定量反演，主要包括声发射P波到时/幅值自动拾取、3D声发射源定位及位移不连续模型震源机制反演3个子程序，代码流程如图 8所示.
|
图 8 位移不连续模型声发射矩张量定量反演代码流程 Fig. 8 Flowchart of the AE source moment tensor analysis program based on the DD model |
本文对定位误差在2.5 mm以内的1149个声发射事件进行震源机制反演，其中692个源事件矩张量反演误差(传感器处的法向测量位移与矩张量计算法向位移差绝对值的加权平均)在20%以内.这692个事件震源反演误差平均值在10%以内，且矩张量特征值均满足约束条件(式(15))，这表明求解结果为仅考虑位移不连续模型的震源机制解，其中部分典型结果如表 3所示.
|
|
表 3 部分典型震源事件矩张量反演结果 Table 3 The moment tensors for some typical AE sources |
根据震源机制反演结果绘制微裂纹模式角α及累计体积的随应变(应力)的变化情况如图 9所示，图中法向体积为位移不连续向量b在裂纹法向的分量bnΔA，而切向体积为位移不连续向量b在裂纹切向的分量btΔA.在整个加载过程中几乎所有的声发射震源机制为混合型，即微裂纹既含有张开/闭合成分也含有剪切滑移成分.为了进一步分析微裂纹的属性，本文定义：当0°≤α < 30°时，微裂纹张开为主；当30°≤α < 60°时，微裂纹为张开和剪切混合型；当60°≤α < 120°时，微裂纹剪切滑移为主；当120°≤α < 150°时，微裂纹为闭合和剪切混合型；当150°≤α < 180°时，微裂纹闭合为主.从图 9a中可以看出，90%以上声发射源由微裂纹剪切滑移产生，其余的为混合型微裂纹，整个过程中没有出现以微裂纹张开/闭合为主的声发射事件，即没有出现微裂纹模式角α介于0~30°、150~180°的事件.
|
图 9 震源参数随应力的变化规律 (a)微裂纹模式角α; (b)微裂纹累计体积. Fig. 9 Variation of AE source parameters with stress based on DD model (a) crack mode angle α; and (b) cumulative microcrack volume. |
在加载初期，微裂纹活动不活跃，声发射事件率低(2-3事件/s)，试样原始裂隙压密闭合震源机制占优(α>90)，微裂纹累积法向体积为负且逐渐减小；随着加载进行，声发射事件率逐渐升高(10-20事件/s)，新微裂纹开始生成并聚集在预制裂纹尖端塑性区内，裂纹模式角α逐渐减小并小于90°，累积切向和法向裂纹体积均出现了较大增长；当应力增大到65.01 MPa(应变0.866%)，裂纹起裂中间岩桥贯通伴随4.12 MPa的局部应力降，声发射事件率突增至约每秒70个事件，出现了10-3mm3量级的大尺度声发射震源事件，累积切向和法向裂纹体积也都呈现出阶梯状抬升，裂纹模式角α降到了最低值；此后，微裂纹活动更加活跃，声发射事件率大幅度增大(30-60事件/s)，绝大多数的裂纹模式角α介于60°~90°间(均值约为80°)，累积切向裂纹体积曲线越来越陡，切向裂纹体积增长速率远大于法向裂纹体积，这些均表明了裂纹起裂后震源事件以微裂纹剪切滑移为主；峰值应力附近，试样处于自组织临界失稳态，大量微裂纹汇合、贯通、成核，伴生着多次局部应力降能量释放，声发射事件率达到最大值(超过100事件/s)，再次出现大体积裂纹震源事件群，震源类型几乎保持不变(α在80°附近波动)，而裂纹累积切向体积分量约为法向体积分量的6倍，同样表明了微裂纹以剪切类型为主；峰值过后，宏观主破裂在瞬间形成(见图 5f)，应力跌落到最低值，由于2个传感器脱落，峰后震源机制解效果不太理想.
4 讨论 4.1 微裂纹空间取向及运动方向确定在3.1节中已经指出，由于矩张量的对称性，微裂纹法向n和位移不连续向量b不能唯一确定(b和n可以互换).在微观尺度上，微裂纹取向很大程度上会受试样内部晶粒分布非均质影响；但在宏观统计尺度上，仍可以认为微观裂纹与宏观主裂纹取向相一致.因此可以根据试样宏观裂纹取向信息来唯一确定b和n，即选取微裂纹法向n使得微裂纹取向与试样宏观裂纹取向尽可能的接近.
按照试样最终破坏模式(图 5f)及宏观主裂纹的取向将试样分为两个区域：(1)区域1(图 10a中黑色宏观主裂纹)裂纹与垂直轴夹角为5.5°左右，该区域裂纹法向为(0.995, 0.096, 0)；(2)区域2(图 10a中红色宏观主裂纹)裂纹与垂直轴夹角为10.5°左右，该区域裂纹法向为(0.983, -0.182, 0).根据声发射定位和震源矩张量反演结果，分别选取与两个区域主裂纹位置接近的声发射事件(区域1中145个事件，区域2中225个事件).需要指出的是，由于其它声发射震源处没有出现明显的宏观主裂纹，无法唯一确定b和n，因此本文暂未对这些震源事件进行分析.根据b和n反演结果绘制微裂纹的空间取向及其运动方向(位移不连续向量b)，结果如图 10b所示.分别统计区域1和区域2中微裂纹的取向及其运动方向与宏观主裂纹的夹角，结果如图 10c和图 10d所示.微裂纹取向与宏观主裂纹取向偏差较大，其夹角的平均值为40.9°(区域1和区域2中两者间夹角平均值分别为43.5°和39.3°).但微裂纹运动方向与宏观主裂纹取向非常接近，其夹角的平均值为17.7°，其中区域1和区域2中微裂纹运动方向与宏观主裂纹夹角平均值分别为18.2°和17.4°.
|
图 10 微裂纹空间取向及其运动方向 (a)破裂模式图; (b)微裂纹及其运动取向; (c)微裂纹与宏观裂纹夹角分布; (d) b向量与宏观裂纹夹角分布. Fig. 10 Orientations of the microcracks and corresponding displacement discontinuity vectors b (a) Failure mode; (b) Microcrack orientations and displacement vector direction; (c) Angle distributions between microcrack orientation and failure plane; (d) Angle distributions between displacement vector b and failure plane. |
在3.5节中，将微裂纹体积沿裂纹法向和切向分解为法向体积分量bnΔA和切向体积分量btΔA，得到了震源以切向滑移为主机制的结论；由于微裂纹法向的任意性，该分解无法反映出全局坐标系下裂纹扩展的断裂力学行为(Kao et al., 2011).在唯一确定微裂纹空间取向及运动方向后，我们可将微裂纹体积在全局坐标系下(坐标系见表 1)分别沿X、Y、Z轴进行分解，则X、Y、Z方向上微裂纹体积分别为Volx=|b|ΔAcosα、Voly=|b|ΔAcosβ和Volz=|b|ΔAcosγ, 其中(cosα, cosβ, cosγ)为位移不连续向量b的方向余弦.在砂岩试样受载过程中，区域1和区域2中的370个震源事件在X、Y、Z方向上累积微裂纹体积分量随应力/应变变化如图 11所示.
|
图 11 全局坐标下累积裂纹体积分量变化规律 Fig. 11 Variation of cumulative microcracks volume components under the global coordinate |
在全局坐标系下，X轴方向微裂纹体积分量Volx对着断裂力学中的Ⅰ型张开或闭合行为，在宏观主裂纹形成及其扩展过程中，Volx呈逐渐增大趋势，在峰值应力处累积Volx分别为累积Voly和Volz的2倍和10倍以上，这表明了受载砂岩试样裂纹扩展中微裂纹事件以Ⅰ型张开机制为主；Y轴方向微裂纹体积分量Voly对着断裂力学中的平面内剪切(Ⅱ型滑开)行为，随着加载的进行，累积Voly逐渐减小，表明Voly主要沿Y轴的负方向，这与试样的全局变形量方向一致，这也是有些学者采用全局变形量来唯一确定b和n的原因(Carvalho, 1999; Chang and Lee, 2004)；而Z轴方向微裂纹体积分量Volz对着断裂力学中的平面外剪切(Ⅲ型撕开)行为，累积Volz在三者中最小，随加载过程在零轴附近波动.虽然试样厚度较小可近似认为该试验为二维力学条件问题，但由于试样内部晶粒分布的非均质性，产生了Ⅰ-Ⅲ型混合加载力学条件(Hull, 1994)，进而造成了朝平面外方向的微裂纹剪切滑移行为.
4.3 模型对比验证
为了验证位移不连续模型的正确性，本文采用震源极性法进行对比验证(Graham et al., 2010).即对于传感器阵列中的任意一个声发射源事件，其震源属性可以通过不同位置处P波震源极性值pol=
|
图 12 震源极性法与位移不连续模型对比 Fig. 12 Variations and comparisons of microcracks type for polarity method and DD model |
试样受载中的局部应力降Δσ也可以通过声发射震源震级来评估，即Δσ=7M0/16/r3，其中M0为震源标量震级(M0= 
(1) 在AIC准则到时拾取基础上，增添了噪声水平统计、到时点验证等功能，提高了声发射定位及震源机制定量反演输入数据的质量；且在背景噪声较大时，该法仍能得到较准确P波到时.采用单纯形最小绝对偏差稳健定位算法对声发射源进行定位，通过定义误差参数来对定位精度定量评价；定位精度在2.5 mm内的声发射事件密集区与试样宏观主裂纹相吻合，故被选为震源机制反演的输入文件.
(2) 针对低孔隙硬岩微破裂震源机制反演问题，建立了仅涉及裂纹张开/闭合和剪切滑移的位移不连续震源模型，即矩张量特征值与试样泊松比之间必须满足M(2)=ν(M(1)+M(3)).通过提出的幅值标量敏感参数获得了传感器处的绝对法向位移值，并采用拉格朗日乘子和Levenberg-Marquardt迭代法求解位移不连续模型下微裂纹的声发射震源矩张量.在上述研究基础上，设计并编写了包含岩石声发射波形检测、3D震源定位和位移不连续模型震源机制反演的岩石裂纹扩展声发射定量反演程序，实现了基于声发射监测的岩石裂纹扩展定量分析.
(3) 含预制裂纹砂岩试样受载裂纹扩展过程中90%以上的微裂纹事件模式角α介于60°~120°，而裂纹累积切向体积分量btΔA约为法向体积分量bnΔA的6倍，表明加载过程中微裂纹扩展以剪切滑移机制占优.选取微裂纹法向n使微裂纹与宏观裂纹取向尽可能的接近来唯一确定微裂纹法向n和运动向量b；微裂纹取向与宏观主裂纹之间夹角平均值为40.9°，而微裂纹运动方向与宏观主裂纹取向非常接近，两者间的夹角平均值为17.7°.在全局坐标系下微裂纹体积分解结果表明，试样裂纹扩展中微裂纹事件X方向上以Ⅰ型张开机制为主，Y方向上代表的Ⅱ型平面内剪切滑移与试样全局变形方向相一致，由于试样内部晶粒分布的非均质性造成了少量微裂纹沿Z轴方向上的Ⅲ型平面外剪切滑移行为.
(4) 对比加载过程中由位移不连续模型获得的裂纹模式角α值与震源极性法的pol值，两者呈一致的变化趋势；根据声发射震源震级及尺度评估试样受载过程中的局部应力降，结果与应力-应变曲线观测到的数据相吻合.以上两点均表明位移不连续模型在砂岩裂纹扩展震源机制定量反演中的适用性.
致谢 感谢美国明尼苏达大学-双城校区Labuz J F教授及美孚公司Kao C S博士在位移不连续模型及震源机制反演理论上的指导和帮助.感谢审稿专家和编辑对论文提出的宝贵修改建议.
Aki K, Richards P G. 1980. Quantitative Seismology: Theory and Methods. Vol. 1. San Francisco: W. H. Freeman and Company.
|
Allen R. 1982. Automatic phase pickers:their present use and future prospects. Bulletin of the Seismological Society of America, 72(6B). |
Bai T Y, Wu S C, Wang J J, et al. 2016. Methods of P-onset picking of acoustic emission compression waves and optimized improvement. Chinese Journal of Rock Mechanics and Engineering (in Chinese), 35(9): 1754-1766. DOI:10.13722/j.cnki.jrme.2015.1565 |
Carpinteri A, Xu J, Lacidogna G, et al. 2012. Reliable onset time determination and source location of acoustic emissions in concrete structures. Cement and Concrete Composites, 34(4): 529-537. DOI:10.1016/j.cemconcomp.2011.11.013 |
Carvalho F C S. 1999. Characterizing brittle failure through quantitative acoustic emission[Ph. D. thesis]. Minnesota: University of Minnesota.
|
Carvalho F C S, Labuz J F. 2002. Moment tensors of acoustic emissions in shear faulting under plane-strain compression. Tectonophysics, 356(1-3): 199-211. DOI:10.1016/S0040-1951(02)00385-2 |
Cerrillo C, Jiménez A, Rufo M, et al. 2014. New contributions to granite characterization by ultrasonic testing. Ultrasonics, 54(1): 156-167. DOI:10.1016/j.ultras.2013.06.006 |
Chang S H, Lee C I. 2004. Estimation of cracking and damage mechanisms in rock under triaxial compression by moment tensor analysis of acoustic emission. International Journal of Rock Mechanics and Mining Sciences, 41(7): 1069-1086. DOI:10.1016/j.ijrmms.2004.04.006 |
Feignier B, Young R P. 1992. Moment tensor inversion of induced microseisnmic events:Evidence of non-shear failures in the -4 < M < -2 moment magnitude range. Geophysical Research Letters, 19(14): 1503-1506. DOI:10.1029/92GL01130 |
Ge M C. 2003. Analysis of source location algorithms:Part Ⅱ.Iterative methods. Journal of Acoustic Emission, 21: 29-51. |
Graham C C, Stanchits S, Main I G, et al. 2010. Comparison of polarity and moment tensor inversion methods for source analysis of acoustic emission data. International Journal of Rock Mechanics and Mining Sciences, 47(1): 161-169. DOI:10.1016/j.ijrmms.2009.05.002 |
Hampton J, Gutierrez M, Matzar L, et al. 2018. Acoustic emission characterization of microcracking in laboratory-scale hydraulic fracturing tests. Journal of Rock Mechanics and Geotechnical Engineering, 10(5): 805-817. DOI:10.1016/j.jrmge.2018.03.007 |
Hull D. 1994. The effect of mixed mode I/Ⅲ on crack evolution in brittle solids. International Journal of Fracture, 70(1): 59-79. DOI:10.1007/BF00018136 |
Ishida T, Labuz J F, Manthei G, et al. 2017. ISRM suggested method for laboratory acoustic emission monitoring. Rock Mechanics and Rock Engineering, 50(3): 665-674. DOI:10.1007/s00603-016-1165-z |
Kao C S, Carvalho F C S, Labuz J F. 2011. Micromechanisms of fracture from acoustic emission. International Journal of Rock Mechanics and Mining Sciences, 48(4): 666-673. DOI:10.1016/j.ijrmms.2011.04.001 |
Kurz J H, Grosse C U, Reinhardt H W. 2005. Strategies for reliable automatic onset time picking of acoustic emissions and of ultrasound signals in concrete. Ultrasonics, 43(7): 538-546. DOI:10.1016/j.ultras.2004.12.005 |
Lennartz-Sassinek S, Main I G, Zaiser M, et al. 2014. Acceleration and localization of subcritical crack growth in a natural composite material. Physical Review E-Statistical, Nonlinear, and Soft Matter Physics, 90(5-1): 052401. DOI:10.1103/PhysRevE.90.052401 |
Li N, Wang E Y, Sun Z Y, et al. 2014. Simplex microseismic source location method based on L1 norm statistical standard. Journal of China Coal Society (in Chinese), 39(12): 2431-2438. DOI:10.13225/j.cnki.jccs.2013.1855 |
Lin Q, Labuz J F. 2013. Fracture of sandstone characterized by digital image correlation. International Journal of Rock Mechanics and Mining Sciences, 60: 235-245. DOI:10.1016/j.ijrmms.2012.12.043 |
Liu J P, Li Y H, Xu S D, et al. 2015. Cracking mechanisms in granite rocks subjected to uniaxial compression by moment tensor analysis of acoustic emission. Theoretical and Applied Fracture Mechanics, 75: 151-159. DOI:10.1016/j.tafmec.2014.12.006 |
Liu P X, Liu L Q, Huang Y M, et al. 2009. Robust arithmetic for acoustic emission location. Chinese Journal of Rock Mechanics and Engineering (in Chinese), 28(S1): 2760-2764. |
Liu P X, Chen S T, Guo Y S, et al. 2014. Moment tensor inversion of acoustic emission. Chinese Journal of Geophysics (in Chinese), 57(3): 858-8. DOI:10.6038/cjg20140315 |
Liu Q S, Liu Q, Pan Y C, et al. 2018. Microcracking mechanism analysis of rock failure in diametral compression tests. Journal of Materials in Civil Engineering, 30(6): 04018082. DOI:10.1061/(asce)mt.1943-5533.0002251 |
Maji A, Shah S P. 1988. Process zone and acoustic-emission measurements in concrete. Experimental Mechanics, 28(1): 27-33. DOI:10.1007/BF02328992 |
Ming H J, Feng X T, Zhang C Q, et al. 2013. Moment tensor analysis of attitude characterization of hard rock newborn fracture surface based on microseismic informations. Rock and Soil Mechanics (in Chinese), 34(6): 1716-1722. |
Mooney H M. 1974. Some numerical solutions for lamb's problem. Bulletin of the Seismological Society of America, 64(2): 473-491. |
Moradian Z, Einstein H H, Ballivy G. 2016. Detection of cracking levels in brittle rocks by parametric analysis of the acoustic emission signals. Rock Mechanics and Rock Engineering, 49(3): 785-800. DOI:10.1007/s00603-015-0775-1 |
Nanjo K Z. 2017. A fiber-bundle model for the continuum deformation of brittle material. International Journal of Fracture, 204(2): 225-237. DOI:10.1007/s10704-016-0175-x |
Othsu M. 1991. Simplified moment tensor analysis and unified decomposition of acoustic emission source:Application to in situ hydrofracturing test. Journal of Geophysical Research:Solid Earth, 96(B4): 6211-6221. DOI:10.1029/90JB02689 |
Othsu M. 1995. Acoustic emission theory for moment tensor analysis. Research in Nondestructive Evaluation, 6(3): 169-184. DOI:10.1007/BF01606380 |
Petružálek M, Jechumtálová Z, Kolář P, et al. 2018. Acoustic emission in a laboratory:mechanism of microearthquakes using alternative source models. Journal of Geophysical Research:Solid Earth, 123(6): 4965-4982. DOI:10.1029/2017JB015393 |
Rice J R. 1980. Elastic wave emission from damage processes. Journal of Nondestructive Evaluation, 1(4): 215-224. DOI:10.1007/BF00571803 |
Scruby C B, Baldwin G R, Stacey K A. 1985. Characterisation of fatigue crack extension by quantitative acoustic emission. International Journal of Fracture, 28(4): 201-222. DOI:10.1007/BF00035216 |
Shah K R, Labuz J F. 1995. Damage mechanisms in stressed rock from acoustic emission. Journal of Geophysical Research:Solid Earth, 100(B8): 15527-15539. DOI:10.1029/95jb01236 |
Silver P G, Jordan T H. 1982. Optimal estimation of scalar seismic moment. Geophysical Journal of the Royal Astronomical Society, 70(3): 755-787. DOI:10.1111/j.1365-246X.1982.tb05982.x |
Wang X R, Liu X F, Wang E Y, et al. 2017. Experimental research of the AE responses and fracture evolution characteristics for sand-paraffin similar material. Construction and Building Materials, 132: 446-56. DOI:10.1016/j.conbuildmat.2016.12.028 |
Wang X R, Wang E Y, Liu X F, et al. 2018. Macro-crack propagation process and corresponding AE behaviors of fractured sandstone under different loading rates. Chinese Journal of Rock Mechanics and Engineering (in Chinese), 37(6): 1446-1458. DOI:10.13722/j.cnki.jrme.2017.1672 |
Wang X R, Wang E Y, Liu X F. 2019. Damage characterization of concrete under multi-step loading by integrated ultrasonic and acoustic emission techniques. Construction and Building Materials, 221: 678-690. DOI:10.1016/j.conbuildmat.2019.06.105 |
Wong L N Y, Xiong Q Q. 2018. A method for multiscale interpretation of fracture processes in carrara marble specimen containing a single flaw under uniaxial compression. Journal of Geophysical Research:Solid Earth, 123(8): 6459-6490. DOI:10.1029/2018JB015447 |
Wu S C, Huang X Q, Chen F, et al. Moment tensor inversion of rock failure and its application. Rock and Soil Mechanics (in Chinese), 37(S1): 1-18, doi: 10.16285/j.rsm.2016.S1.001.
|
白添羊, 吴顺川, 王进进, 等. 2016. 岩石破裂声发射压缩波到时拾取方法及其优化改进研究. 岩石力学与工程学报, 35(9): 1754-1766. DOI:10.13722/j.cnki.jrme.2015.1565 |
李楠, 王恩元, 孙珍玉, 等. 2014. 基于L1范数统计的单纯形微震震源定位方法. 煤炭学报, 39(12): 2431-2438. DOI:10.13225/j.cnki.jccs.2013.1855 |
刘培洵, 陈顺云, 郭彦双, 等. 2014. 声发射矩张量反演. 地球物理学报, 57(3): 858-866. DOI:10.6038/cjg20140315 |
刘培洵, 刘力强, 黄元敏, 等. 2009. 声发射定位的稳健算法. 岩石力学与工程学报, 28(S1): 2760-2764. |
明华军, 冯夏庭, 张传庆, 等. 2013. 基于微震信息的硬岩新生破裂面方位特征矩张量分析. 岩土力学, 34(6): 1716-1722. |
王笑然, 王恩元, 刘晓斐, 等. 2018. 裂隙砂岩裂纹扩展声发射响应及速率效应研究. 岩石力学与工程学报, 37(6): 1446-1458. DOI:10.13722/j.cnki.jrme.2017.1672 |
吴顺川, 黄小庆, 陈钒, 等. 2016. 岩体破裂矩张量反演方法及其应用. 岩土力学, 37(S1): 1-18. DOI:10.16285/j.rsm.2016.S1.001 |