波形分析是指以声发射信号的时域波形为对象，选用一定的信号处理方法来获取声发射源信息。声发射波形信号中蕴含着大量声发射源的信息。因此，近年来，波形信号分析成为了声发射技术研究的新热点。声发射波形分析方法主要有基于傅里叶变换的频谱分析和时频分析两种方法。在频谱分析方面最突出的研究成果，是在室内实验[13-14]与现场监测[15-16]中发现了不同尺度岩石破坏时出现的声发射信号由高频向低频移动的“频移”现象。傅里叶变换的不足之处是信号的“整体变换”，即不能揭示频率分量随时间变化的情况。因此，只适用于平稳信号的分析[17]。时频分析(time-frequencyanalysis)的基本思想是设计时间和频率的联合函数，用于描述信号在不同时刻的谱分量[18]。时频分析将一维时域信号变换为二维时频图像，能够展现出更为丰富的信号内涵。时频分析包括短时傅里叶变换和小波变换两类线性时频变换，和以维格纳变换为代表的双线性时频变换[18]。小波分析是声发信号多尺度分析的有力工具[19-20]。然而，小波变换给出的是时间–尺度(或分辨率)之间的关系[19]，对于要求频率信息的场合，时频变换具有不可替代的作用。时频变换在岩石破坏前兆信号分析方面有诸多应用：例如，对岩石试件多轴加载的声发射信号进行分析，发现低应力水平时声发射信号为高频小幅值、岩爆或地震前兆为低频大幅值。
与数字信号处理的其他领域相比，岩石破坏声发射信号时频分析的研究还很不充分，现有的岩石破坏频域前兆包含的信息仅有时间、频率和幅度，缺少动力学信息的描述；二次型时频分布的“交叉项”识别以及频率的精确定位缺少综合性的方法，因而与实际应用的要求还有不小的距离，其主要问题是：(1)大数据：即声发射原始波形的采样为海量数据，必须分割后才能进行处理；对大数据的分割的研究少有报道。(2)对声发射信号滤波及时频分析算法的研究不够充分。(3)对时频图像的进一步处理与分析的研究，尚不能达到细致描述岩石破坏动态过程的要求。本文以花岗岩单轴压缩声发射数据为实例，开展大数据样本的分割理论、时频分析窗函数滤波、时频分析算法、和以时–频联合分布、谱分量分布模式表征动态过程为内涵的，多维度瞬时频率前兆信息的研究，以期为岩石破坏声发射信号前兆信息的识别提供原理上的支持。
实验系统及原理
岩石试样为山东莱州花岗岩，呈浅肉红色，粗粒块状结构，肉眼观察试件致密、无裂纹(见图1(b))。花岗岩是典型的硬岩，由于工程开挖，在高应力作用下，会产生突然猛烈的破坏。将花岗岩块首先用岩石切割机进行切割，再用双面磨石机磨平，制成规格为150mm×60mm×30mm(高×长×宽)的长方形试样，试样两侧安装有2个声发射传感器，直接固定在试样表面(见图1(a))。矿物成分的X射线衍射(X-raydiffraction，XRD)分析表明，花岗岩中晶体矿物主要为石英(27.0%)、钾长石(37.0%)、钭长石(31.0%)和黏土矿物(5.0%)。扫描电镜(scanningelectronicmicroscope，SEM)观察分析表明：花岗岩是非均匀、多孔性介质；石英、长石、云母接触边界明显可见；石英、长石中存在着大量的微米级孔洞。岩石的非均匀性源自于其中存在着不同硬度的成分、存在着微观孔洞与不同成分的边界。当岩石加载时，这些微观缺陷会扩展、传播、最终汇合成宏观裂纹并导致岩石的完全破坏。实验是在中国矿业大学(北京)煤炭资源与安全开采国家重点实验室进行，实验系统如图2所示。花岗岩试件的单轴加载实验是在EHF–UG500全数字液压伺服试验机进行，最大静态载荷为759kN，荷载精度小于0.5%，采用轴向位移控制，加载速率为0.002mm/s，自动记录实验过程中的力和位移数据。在对花岗岩样的加载过程中，同步进行声发射探测与红外探测。本文主要研究声发射信号的处理。声发射探测采用中国矿业大学(北京)深部岩土与地下工程国家重点实验室的PXWAE声发射监测系统。系统由传感器、放大器、信号采集处理和记录与显示(采集卡和微机)组成(见图3[2])，最高采样速率为20MHz，采样精度12位。采集卡与计算机的数据通讯采用PCI总线，可实现速率达132MB/s的波形信号传输。传感器接收从试件内发出的波动信号，使声发射传感器压电陶瓷变形，电压发生变化。这种变化的电压信号通过前置放大器放大，再经过信号调制变成稳定的声发射信号存储至微机中并实时显示波形等参数。实验采用的声发射传感器的共振频率为150kHz，前置放大器40dB，增益设置为10倍。实验中，系统的采样率设置为1MHz。
声发射信号分割
1.信号分割理论
图4给出了一个拉长了声发射信号随时间t的发展过程(基本波形)，代表了由一个特定尺度的破裂事件引起的应力波，对应于岩石内部热力学与动力学状态的变化过程。对基本波形的描述，除了在节1中介绍的统计参数外，还可以将其看作一个基本时间序列，用幅度A(基本波形最大值)、基本周期T(基本波形中2个相邻峰值间的时间)、基本频率f(基本周期的倒数f=1/T)、能量(定义为波形包络线下的面积)、以及持续时间D(与波形的能量成正比)等来描述其动态特征。声发射信号的频域特性与震源类型、震源与传感器之间的距离有关，单个声发射基本波形可能含有不同的频率成分。声发射信号的频谱分布取决于2个因素：(1)震源信号的频谱；(2)传播到传感器的过程中所经历的衰减(或改变)。信号衰减对信号的频谱分布起着重要作用。一般的规律是：衰减随频率的增加而增加。因此，在较长的距离，只有低频分量的声发射信号可观测到。并且，如震源信号不含显著的低频信号时，则存在着一个临界距离，在此距离之外，此震源的声发射信号不能被检测到。图5给出了声发射信号相关研究的频率分布范围。岩石中产生声发射信号的源有很多，主要有位错、胶结物破裂、相变、孪生/滑移变形等；断裂扩展的不同阶段，即微裂纹起裂、扩展与汇合[1]。因此，声发射信号在不同频带上的分量可以反映岩石破坏的多尺度信息。相对于时域声发射统计参数，其基本周期和基本频率具有统计性质上的确定性。随着计算机技术的进步，为了探测到微裂纹起裂与扩展过程，现今岩石力学实验室的声发射探测大都采用约1MHz以上的采样频率。通常一个岩石力学实验，如岩爆实验需要10～30min的时间，所得到的声发射数据的点数(样本数)为海量级，形成了所谓的“大样本数据”。常见的大样本数据来源于语音学及各种地球物理现象的表达。目前常用按固定时间进行分段(时窗分割)后再进行傅里叶变换处理[24]。时窗长度的选取要比原数据长度充分小，以便能够描述快速变化的信息。同时，时窗又要足够长，使其中包含的信息具有因果关系。图6给出了大样本数据的时域和频域分割原理。其中x[n]为N点时间序列的大数据，n为离散时间变量。设采样频率为fs(Hz)，则大数据的基本频率为0maxff/N(1)式中：maxsff为大数据的采样频率(Hz)，也是最高频率，其倒数ssT1/f为采样周期(s)，即图4中的基本周期T。设以时窗宽度为wN对大样本进行分割，wN也表示时窗数据点数。则时窗宽度(时间分辨率)与频窗宽度wf(频率分辨率)的关系为wmaxwff/N(2)设信号x[n]的离散傅里叶变换为F()，为离散圆频率，与时域分割相对应的频域分割如图6(b)所示。窗宽越小，相应域内的分辨率就越高。根据Heisenberg测不准原理(wwNf≥(1/(4))，时间分辨率与频率分辨率不能同时任意小，它们的乘积受到一定的限制，要提高时间分辨率就要降低频率分辨率。不确定性原理表明：同时有任意小的时宽和任意小频宽的图像是不存在的。因此，对大数据进行时域分割，等于降低了原来的频率分辨率，由式(1)和(2)与图6给出的关系可知，短数据的频率分辨率wf与大数据的基本频率0f的关系为：0wwf(N/N)f，其中wN＞N。即与大数据相比，短数据的分辨率降低了wN/N倍。因此，在进行数据分割时，应当考虑适当增加时间窗口宽度，以提高分割后短数据时频图的频率分辨率。#p#分页标题#e#
2.信号分割实验
在大数据的分割中，时窗wN的选择是关键问题。与本文引言中综述的声发射信号的研究方法相适应[3]：对于参数分析法，由于声发射的参数统计是以图4中的基本波形为对象，因此时窗的长度应与基本波形的长度相同；对于波形分析法，其研究对象是由声发射时域波形表征的因果时间序列，因此，考虑一段时间内的声发射信号更有实际意义。岩石破坏过程经历了裂纹的初生、扩展传播、汇合及至宏观破坏的动态过程[25]。也就是说，岩石的破坏是与应力路径(加载历史)相关的非线性过程。反映声发射信号的波形分析法，其研究对象应当是由一系列基本波形组成的具有一定因果关系的时间序列，能够反映材料在一定应力水平下裂纹扩展的动力学过程。因此，当采用波形分析法时，声发射大数据分割的时窗长度的选择应当根据加载方式(动载荷、静载荷)、岩样尺度、岩样物理性质等因素来选择。如动载荷、应选用较短的时窗；静载荷情况下应选用较长的时窗。另外，对大数据分割后的长度为wN的短数据，采用短时傅里叶变换进行分析时，变换本身的窗函数h(n)的长度应远小于wN，即能够反映短数据的时频局部化性质。事实上，由于问题的复杂性，短数据长度wN的选择存在多种可能性。除去上述原则外，还要由大量的声发射信号数据处理实验来决定。在深部岩土力学与地下工程国家重点实验室，对于不同地区的、同样尺寸的花岗岩试件进行了百例以上的应变岩爆实验与几十例单轴压缩实验[1，21]。
其声发射采样频率均为为1024kHz，信号长度在109～1010数量级。为了开展时频分析研究，将声发射大数据用不同的时窗长进行分割，最后选定8ms为时窗长度。本文以46#花岗岩单轴压缩声发射信号为研究对象，采样频率为1024kHz，信号总长为264192000点，分割时窗长为8ms。分割后声发射短数据的3种主要模式如图7所示。图7(a)为由一个基本波形组成的单震型破坏事件，多发生在加载初期；图7(b)为由2个基本波形组成的所谓前震–主震型破坏事件[14]；图7(c)为由2个双前震–主震型组成的多震型破坏事件。实际声发射短数据的模式比较复杂，存在连续的基本波形组成的事件。限于条件，目前笔者只能选用固定的时窗长度，用硬件进行分割。值得一提的是，大数据分割的时窗长度应当是自适应的，即能够容纳一个或多个完整的基本波形。本文选用8ms时窗长的原因之一是，分割后的短数据中不完整的波形最少。尽管如此，在分割后的短数据中，仍有图7(c)中的最后一个波形的一小部分没有完全显现的情况。然而，时频算法是对数据分段进行变换。作为补救措施，对于应用短时傅里叶变换类的时频变换，应当使时频变换中分析窗的长度远小于短数据长度，来反映信号的局部性质，以使丢掉波形中的信息损失最小；当采用双线性时频法时，由于信号的自相关长度远小于时窗长，尾部一小部分的丢失的信息并非短数据波形的主要成分。
时频分析基本理论
典型的线性时频变换(或线性时频表示)是短时傅里叶变换(short-timeFouriertransform，STFT)和Gabor变换(Gabor可以看作是最优短时傅里叶变换)。线性时频表示是由傅里叶变换演化而来的，满足线性叠加性。1946年Gabor提出的短时傅里叶变换实质上是加窗的傅里叶变换[26]，即为了达到时域上的局部化，在信号傅里叶变换前乘一个时间有限的窗函数h(t)，并假定非平稳信号在分析窗的短时间隔内是平稳的，通过窗函数h(t)在时间轴上的移动，对信号进行逐段分析得到信号的一组局部“频谱”。信号x(t)的短时傅里叶变换定义[26]为2π()()()edjfSTFTtfxht，(3)式中：h(t)为分析窗函数。式(3)表明，信号x(t)在时间t处的短时傅里叶变换就是信号乘上一个以t为中心的“分析窗”h(t)后所作的傅里叶变换。x(t)乘以分析窗函数h(t)等价于取出信号在分析时间点t附近的一个切片。对于给定时间t，STFT(t，f)可以看作是该时刻的频谱。
特别是，当窗函数取h(t)1时，则短时傅里叶变换就退化为传统的傅里叶变换。要得到最优的局部化性能，时频分析中窗函数的宽度应根据信号特点进行调整，即正弦类信号用大窗宽，脉冲型信号用小窗宽。而STFT的窗宽是固定的，不能进行自适应调整；但其最大优点是其基本算法即是傅里叶变换，易于解释其物理意义，不存在双线性时频表示的交叉项问题。双线性时频表示也称作二次型时频表示，其代表是Cohen类时频分布。它反映信号能量的时频分布。二次型时频表示不满足线性叠加性。设信号为12x(t)ax(t)bx(t)，记x(t)，1x(t)，2x(t)的线性时频变换分别为P(t，f)，1P(t，f)，2P(t，f)，则有如下关系[26]：2212P(t，f)|a|P(t，f)|b|P(t，f)122R[abP(t，f)](4)式(4)中等号右边最后一项称为交叉项(干扰项)，是二次型时频表示的固有属性。信号x(t)的Cohen类时频分布的一般表达式[26]为（略）式中：(，v)为核函数，x为x的共轭函数。Cohen类中最重要的成员是1932年Wigner提出了维格纳分布，最初用于量子力学的研究。1948年，Ville将其引入信号分析领域，称为Wigner-Ville分布(Wigner-Villedistribution，WVD)[27]。在式(3)中取核函数(，v)1，信号x(t)的WVD定义[26]为（略）由定义可知，WVD(t，f)为实值函数，在时频平面中分别是信号能量和功率谱分布，为信号自相关的傅里叶变换，不存在STFT中的固定窗口问题，具有很好的时频聚集性，是分析非平稳时变信号的有力工具。然而，对于多分量信号，根据卷积定理，WVD(t，f)会出现交叉项，如式(4)所示。交叉项会产生“虚假信号”，这是WVD(t，f)应用中的主要缺陷。交叉项通常是振荡的，而且幅度可以达到自相关项的2倍，造成信号特征模糊不清。与短时傅里变换不同。应用Cohen类(二次型类)时频分布应用的主要问题是如何抑制交叉项[23]。人们进行了许多关于Cohen类时频分布交叉项抑制或消除的算法的研究[26-27]：交叉项与时频分布的有限支撑特性密切相关，而交叉项的抑制又主要通过核函数的设计来实现。典型的加核函数的WVD有伪维格纳分布(Pseudo-Wigner-Ville，PWVD)和平滑伪维格纳分布(smoothedpseudoWigner-Villedistribution，SPWVD)等等。SPWVD是近年来受到重视的一种抑制交叉项的算法，其定义[26]为（略）式中：g(s)为平滑核函数，h()为窗函数[28]。
从一个连续信号中抽取一段，如果不加窗，实际上是默认加了一个矩形窗，用矩形窗函数作为分析窗进行短时傅里叶变换，会使信号的频谱不只是出现在中心频带内，而且在中心频带以外都会有谱分量出现，即出现所谓的“频谱泄漏”现象[29]。由于矩形窗具有旁瓣高、衰减慢的特点，其频谱泄露现象较为严重，虽具有高频率分辩率，但其频谱泄露现象会对时频分析的准确性造成较大的影响。因此，窗函数的设计对于保证时频分析结果的正确性具有重要意义。窗函数设计包括2个问题，一是窗长、二是窗函数滤波性能。窗函数长度涉及分段截取信号，如前所述，这主要根据信号的性质和计算机实验来决定，以能够识识别信号中的最小脉冲为原则，即窗宽应小于等于最短的基本波形的持续时间D。窗函数h(n)一般是低通滤波器。图8给出了典型窗函数的频率响应jH(e)。理想窗函数的选取原则是：(1)窗函数的主瓣应尽量窄，使能量集中在主瓣内，在一定的时窗长度下获得较高的频率分辨率；(2)旁瓣高度尽可能小且随频率增高快速衰减，以减少频谱的泄露失真。关于窗函数的设计，涉及到有限冲击响应滤波器(FIR滤波器)的设计理论，详见V.K.Ingle和J.G.Proakis[29]的研究成果。窗函数的评价指标通常由分贝(dB)形式给出（略）#p#分页标题#e#
计算机数据处理
1.窗函数特性实验
实验利用Matlab8.0时频分析工具箱TFTB2002实现。TFTB2002中提供了许多标准窗函数，如矩形窗(Rectwin)、海明窗(Hamming)、布莱克曼窗(Blackman)、汉宁窗(Hanning)、凯塞窗(Kaiser)、三角窗(Triang)、切比雪夫窗(Chepwin)等等。如前所述，短时傅里叶变换的分析窗长应远小于声发射短数据的长度。本文处理的是8192点声发射短信号，根据声发射基本波形的的最小长度，经大量统计与反复试算，当时窗wN的长度为64点时得到了时间分辨率与频率分辨率的最优妥协效果。海明窗的时域[29]可表示为h(n)0.540.46cos(2πn/N)(0≤n≤N)(9)海明窗的64点时域信号如图9(a)所示，其傅里叶幅度谱由64点DFT(discretefouriertransform)计算，如图9(b)所示。
为了比较，图9(c)和(d)分别给出了64点矩形窗的时间函数与傅里叶谱。根据图8中给出的指标对图9进行分析可知，海明窗的主瓣宽为0.06(rad)，矩形窗的主瓣为0.15(rad)，说明矩形窗有更好的频率分辨率；但海明窗的旁瓣衰减快于矩形窗，说明海明窗具较小的频谱泄漏。一般来说，主瓣即窄、旁瓣小且衰减快的窗函数是不容易得到。因此，为减少频谱泄露现象以获最更真实的谱特征，对文中列出的窗函数进行了一一比较(略去比较结果)，最终选用频谱泄漏小且频率分辨率尚可的64点海明窗作为短时傅里叶变换的最优窗函数。当然，由于不同物理过程时间序列的复杂性，窗函数的选取是与实际问题紧密相关的。
2.时频分析算法实验
应用由式(3)，(6)，(7)定义的STFT，WVD和SPWVD等3种算法进行花岗岩单轴压缩声发射短信号的分析，对频率局部化性能、可识别性进行比较。图10为在0.98c的声发射短信号的Gabor变换(即最优STFT，以下简称为STFT)时频谱图。在较高频率7.2105和8.8105Hz的连续高频分量，以及在3105～4105Hz范围内的连续的低频分量，给出了在此应力水平上的破坏前兆信息。如前所述，STFT是对时间信号进行分段的傅里叶变换，其意义明确。不足之处是频率分辨率较低，即得到的具有峰值的分量分布在一个较宽的频带上，难以确定所谓的“主频”或“中心频率”的准确数值。图11给出了同一个声发射信号的WVD时频谱图。其频率分辨率相对较高，很容易读出主频的数值。然而，如果没有图10中的STFT时频图的先验知识，很难分出哪些是主频、哪些是交叉项形成的分布。图12给出了SPWVD时频谱，其交叉项被显著抑制了，但仍存在交叉项，并且在3105～4105Hz范围内的有用的频率成分也同时被删掉了。
因此，综合上述研究，针对岩石力学实验中的声发射大数据，提出如下声发射信号综合优化算法：(1)首先应用节3中的大数据分割理论对声发射采样数据进行最优分割；(2)再用节4中的窗函数优化理论对窗函数进行优化，用优化了的STFT进行主频带分布的定位；(3)再利用WVD分辨率高的特点，确定主频带的中心频率；(4)对时频分布图进行精细化处理，以便于分析和理解其中蕴含的动态信息，得到岩石破坏的多维度瞬时频率前兆。这种综合优化时频分析方法，可以在准确定位中心频率位置的同时，避免误判，提高分析的可信性。下面考察优化得到的基于64点海明窗的短时傅里叶变换捕捉声发射短信号瞬时频率的性能。应用STFT算法对图7中3个典型的短声发射信号进行变换，其时频分布如图13所示。STFT时频谱图准确反映了图7中各个脉冲频率随时间的变化关系，说明算法具有很好的时频局部化性能。
实验结果分析
花岗岩单轴压缩实验的全应力–变应变曲线如图14所示。点A，B，D为曲线上斜率发生变化的点(拐点)，将曲线分为以下阶段[25]：(1)OA段：为非线性压密阶段；此阶段试件刚度逐渐增加，曲线上弯，煤岩内部的天然缺陷(微裂纹、微空隙等)在外载荷作用下逐渐闭合。(2)AB段：为线弹性变形阶段，点A为微裂纹初生的点；此阶段应力与应变成正比，花岗岩的弹性模量由AB段曲线斜率确定，根据回归分析，其弹性模量为14.2GPa。(3)BC段：为非线性变形阶段，可进一步分为：①裂纹稳定扩展阶段(BB1)：裂纹端部引起局部应力集中或剪切面滑移，促使裂纹稳定扩展，剪切带扩大，应力–应变曲线下弯、刚度减小；②非弹性变形破坏阶段(B1C)：随载荷增加，变形非线性地增长，微裂纹迅速扩展并汇合。点C的应力称为峰值应力或单轴抗压强度(c=114.2MPa)。在点C，微裂纹扩展、汇合成宏观裂纹，岩样破坏。(4)CD段为应变软化阶段，试样达到峰值载荷后，载荷随变形的增长而减小，岩石内大量新微裂纹产生、扩展、汇合导致试件完全破坏，点C1为岩样完全破坏点。
花岗岩的单轴压缩曲线有明显的残余载荷阶段。(5)点D为峰后应力–应变曲线上的拐点，表示试件在较低的应力水平上又重新具有一定的承载能力。研究表明，声发射信号的频率分量与破裂尺度成反比，即高频对应于小尺度的裂纹[30]。在微观尺度下，声发射来源于晶格破坏或错位；在细观尺度下，声发射来源于岩石颗粒转动或边界移动、或沿晶或穿晶破裂；声发射同样可以反映宏观尺度结构的破坏[1，31]。图15给出了应用本文提出的声发射信号综合优化算法，得到的关键点A～D的STFT时频图(限于篇幅，WVD时频图没有给出)。在数据处理过程中，对二维时频图进行精细化处理包括：增加插值层数、提高图像分辨率和伪彩色处理等。分析中，主频(高强度频率的分布频带)由STFT识别，中心频率(主频带中心)由WVD定位。图15(a)给出了对应于应力水平(0.12c)的STFT时频图。可以看到有一个在中心频率6.5105和8.5105Hz分布的高频脉冲型应力波(单震型)，其持续时间为1ms左右。在声发射数据的分析处理中观察到，在花岗岩的压密阶段，8ms采样的声发射信号均为单个脉冲型。表明压密阶段是单个微尺度随机破裂事件。图15(b)为弹性阶段起点A(c0.35)的STFT时频图。有4个接续的短脉冲分布在中心频率6.5105和8.5105Hz。说明此应力水平是微尺度破坏，但破坏是由多个微裂纹的连续破坏形成，为多震型，代表微裂纹的扩展过程。图15(c)为弹性阶段应力水平为(c0.43)的STFT时频图，为前震–主震型，但2个脉冲型分量的持续时间要更长，其幅度更高，同样分布在中心频率为6.5105和8.5105Hz的中高频区附近，代表了弹性阶段微裂纹的持续扩展。图15(d)分别为弹性阶段应力水平c0.50和c0.66的STFT时频图，两图均出现了分布在中心频率3105和4.5105Hz、持续时间短的低强度低频分量，以及在6.5105Hz的较高频率的分量，表明微裂纹的扩展与较大尺度裂纹出现。图15(d)和(e)的时频分布模式，即低强度短持续时间低频分量，可看作是岩石破坏的早期预警信号。在最高频率(106Hz)的分量应视为应力波与机器振动耦合形成的噪声分量。图15(f)(和(g)分别为弹性阶段应力水平c0.78和c0.82(弹性阶段即将结束)的时频图。与图15(d)和(e)相比，在3105和4.5105Hz的低频脉冲分量的持续时间更长、强度更高，表明在裂纹扩展的速率加快；在6.5105Hz附近的较高频率的分量的强度与持续时间均有所降低，表明两图代表的时间段内，裂纹的汇合是主要破裂事件。图15(f)和(g)的时频分布模式，即持续时间较长的高强度低频分量，可看作是岩石的早期破坏前兆信号。图15(h)为弹性阶段结束点B(c0.89)的STFT的时频图。#p#分页标题#e#
在高频带7.2105和8.8105Hz分布着2个长持续时间、较高强度的高频分量，说明大量微裂纹的出现，裂纹扩展速度加快。在低频带3105Hz附近的低强度脉冲型频率分量也变为长持续时间信号，表征了裂纹扩展、成核的动力学过程。图15(f)的时频分布模式应视为破坏前兆。图15(i)为非线性变形阶段开始点B1(c0.98)的STFT时频图。在高频带7.2105和8.8105Hz分布的2个高强度、长持续时间分量，以及在低频带3105Hz分布的高强度、长持续时间的分量，代表了破坏过程的进一步加速、宏观裂纹的出现，是临近破坏的前兆。图15(j)为峰值应力点C(c)的STFT时频图。在中心频率8.8105Hz分布的高强度、长持续时间高频分量仍存在，表示微裂纹的扩展仍在加速；而在低频带3.3105Hz的长持续时间分量的强度达到最大，表明了宏观裂纹的贯通。图15(k)为峰后点C1(c0.95)的STFT时频图。在中心频率6.5105与3105Hz分布的2个分量变成了低强度、短持续时间的应力波。表明峰后的非连续、塑性破坏远不如峰前的破坏剧烈。图15(l)为峰后残余载荷阶段点D(c0.63)的STFT时频图。与峰前的时频图不同，点D仅有在中心频率6.5105Hz上分布的脉冲型分量。表明破碎岩石中仅有随机的小尺度破坏事件，其力学行为发生了根本性的变化。由图4可知，以一定频率传播的声发射波形的幅度与持续时间与波形能量成正比。因此，通过对时频图中的谱分量的幅度与持续时间的分析，可以揭示岩石破坏的机制。精细化处理的时频图不仅给出了时间–频率联合分布，同时也给出了应力波分量不同的分布模式。例如在图15中，低频分量的分布模式有：小幅度短持续时间(图15(d)，(e)，(k))、较高幅度并较长持续时间(图15(f)，(g))、低幅度长持续时间(图15(h))、高幅度长持续时间(图15(i))、高幅度连续分布(图15(j))等5种模式；高频分量的分布模式包括：脉冲型(图15(a)，(l))、前震–主震型(图15(c))、双前震–主震型(图15(b))、较高幅度较长持续时间(图15(d)，(e))、低幅度较长持续时间(图15(f)，(g)，(k))、连续的低幅度双高频分量(图15(h))、连续的高幅度双高频分量(图15(i))、连续的高幅度高频分量和低幅度高频分量(图15(j))等8种分布模式。时频图中谱分量的分布模式及其组合、加上时间–频率联合分布，提供了岩石破坏的多维度瞬时频率前兆信息。
结论
给出了声发射信号波形分析法的大数据优化分割原理、大数据与分割短数据的频率分辨率计算公式。对于波形分析法，声发射时域信号的分割应表达了一个动态过程，分割长度应当表达在一定应力水平上、一个具有动力学因果关系的破坏事件。大数据分割时，还应考虑过短的分割时间窗会降低短数据的频率分辨率、以及数据长度应当适合于个人计算机处理等问题。短时傅里叶变换的分析窗宽度关联于时频局部化性能、分析窗函数性质相关于滤波效果与抑制频谱的泄漏(减少分析结果中的假频率)。分割窗函数宽度应根据信号的性质与计算机实验来决定，以能够识识别信号中的最小脉冲为原则，即窗宽应小于等于最小的基本波形的长度；同时，尽量采用大窗宽，以提高频率分辨率。分析窗函数的优化可由其频率特性来进行，其原则是：(1)窗函数的主瓣应尽量窄，使能量集中在主瓣内，在一定的时窗长度下获得较高的频率分辨率；(2)旁瓣高度尽可能小且随频率增高快速衰减，以减少频谱的泄露失真。
针对岩石力学实验中的声发射大数据，提出如下声发射信号的综合优化算法：(1)首先应用节3中的大数据分割理论对声发射采样数据进行最优分割；(2)再用节4中的窗函数优化理论对窗函数进行优化，用优化了的STFT进行主频带分布的定位；(3)再利用WVD分辨率高的特点，确定主频带的中心频率；(4)对时频分布图进行精细化处理，以便于分析和理解其中蕴含的动态信息，得到岩石破坏的多维度瞬时频率前兆。这种综合优化时频分析方法，可以在准确定位中心频率位置的同时，避免误判，提高分析的可信性。时频分析综合优化算法，对花岗岩单轴压缩声发射短数据进行了分析。时频分布图中的幅度与持续时间给出了谱分量的不同分布模式，可以对裂纹扩展的动力学过程进行详细描述。这些分析结果可以归纳成为以时–频联合分布、谱分量分布模式表征动态过程为内涵的多维度瞬时频率前兆信息。现将本文分析的花岗岩单轴压缩的多维度瞬时频率前兆信息简述如下：(1)中心频率为6.5105和8.5105Hz的脉冲型高频分量为压密的特征。(2)c0.35应力水平上的6.5105和8.5105Hz的多脉冲型高频分量为弹性阶段起点。(3)c0.78水平上3105和4.5105Hz的低频长时间脉冲型高幅度分量，是早期破坏前兆。(4)c0.89水平上7.2105和8.8105Hz的长持续时间高强度高频分量，以及3105Hz低强度长持续时间低频分量，为破坏前兆。(5)c0.98水平上7.2105和8.8105Hz的2个高强度、长持续时间分量，与3105Hz的高强度、长持续时间低频分量，是临近破坏的前兆。(6)峰值应力c的8.8105Hz高强度、长持续时间高频分量和3.3105Hz的长持续时间最高强度的低频分量表明了宏观裂纹的贯通。（本文图略）
本文作者：宫宇新 何满潮 汪政红 尹雨婷 单位：武汉中南民族大学 数学与统计学学院 北京国矿业大学 深部岩土力学与地下工程国家重点实验室
