数字信号处理的几个前沿课题

—PAGE306—数字信号处理的几个前沿课题前面介绍了数字信号处理的基本知识，本章我们将介绍时谱分析、小波变换、地震观测系统仿真与地面运动恢复等几个数字信号前沿课题，以便大家在实际工作中参考。10.1时谱（倒谱）分析时谱分析(Cepstrumanalysis)是一种非线性信号处理技术，它在语言、图像、和噪声处理领域中都有广泛的应用。时谱可分为两类：复时谱和功率时谱。MATLAB信号处理工具箱提供复时谱分析的工具函数。复时谱（Complexcepstrum）的定义为：（10-1）由上式可见，复时谱实际上是序列x(n)的Fourier变换取自然对数，再取Fourier逆变换，得到的复时谱仍然是一个序列。也就是说，复时谱是x(n)从时间域至频率域、频率域至频率域、频率域至时间域的三次变换。MATLAB信号处理工具箱函数cceps用于估计一个序列x的复时谱，调用格式为：xhat=cceps(x)式中，x为输入序列（实序列）；xhat为复时谱（复序列）。MATLAB信号处理工具箱还提供了序列实时（倒）谱的计算程序rceps，调用格式为Y=rceps(x)，其中x为实序列；y为实时谱，执行的操作为：（10-2）由此可知，我们不能从序列x的实时谱重构原始序列，因为实时谱是根据序列Fourier变换的幅值计算的，丢失了相位方面的信息。但如果需要，可采用最小相位模式估计原始序列。由于复时谱从复频谱计算得到，不损失相位信息，因此复时谱是可逆的，实时谱过程是不可逆的。时谱分析技术广泛地应用于语言信号分析、同态滤波技术中。这里举一个说明复时谱在具有回声信号测量中的应用。【例10-1】设原信号是一个45Hz的正弦波，在传播过程中遇到障碍产生回声，回声振幅衰减为原信号的0.5，并与原信号有0.2s的延迟。在某测点测到的信号是原信号和回声信号的叠加。试使用复时谱分析该测点的信号。%Samp10_1t=0:0.01:1.49;%时间信号，采样间隔为0.01ssig=sin(2*pi*45*t);%原始信号echo=0.5*[zeros(1,20),sig(1:length(t)-20)];%在信号前面补加20个零，并使其振幅衰减一半作为回声信号sigecho=sig+echo;c=cceps(sigecho);%求原始信号与回声信号叠加的复时谱subplot(2,1,1),plot(t,sigecho);%绘制原始信号与回声的叠加信号xlabel('时间/s');axis([01.5-11]);gridontitle('时间域信号')subplot(2,1,2),plot(t,c);%绘制其复时谱xlabel('时间/s');axis([01.5-11]);gridontitle('时谱域信号')程序运行结果为图10-1。可以看到，信号的复时谱在t=0.2秒处有一个峰值，这就是回声信号。表示运用复时谱可以明确检测到回声信号。在数字地震仪监测核爆破的工作中，如果采用两次核爆在几乎同一时间、同一地点发生，爆破产生的地震波经过地球内部传到同一地震台上的传播路径几乎一致，并且爆破源产生的波形相差不大，因此在地震台上表现为顺序到达的两个波形类似的波，若后面的波形被淹没，可以采用复时谱的方法进行检测，明确给出是否是两次爆破事件，如果是两次爆破事件，就可以明确给出两次事件的时间间隔。图10-1正弦波与回声信号叠加的波形和时谱形状10.2地震观测系统的仿真和地面运动的恢复由前面的线性系统的讲解可知，一个系统可以用它的系统函数或脉冲响应来表示：（10-3）式中，x(t)为输入信号，对于地震观测系统来讲为地面运动；y(t)为系统的输出，对于地震观测系统来讲为地震记录；h(t)为系统的脉冲响应。在频率域内，根据卷积定理，该式可以表示为：（10-4）式中，为系统的传递函数，和分别为x(t)和y(t)的Fourier变换。现在设想一个频带范围很宽的线性系统，比如宽频带地震仪，其系统函数为;另一个频带较窄的系统，比如短周期地震仪，其系统函数为，对于同样的输入有：，（10-5）式中，为频带较窄的系统记录的频谱；为频带较窄系统的传递函数。比较前面两式可得：（10-6）将上式变换到时间域就得到频带较窄系统的输出y1(t)。也就是说，如果知道了宽频带和窄频带系统的传递函数和，原则上可以从宽频带系统的输出推测出频带较窄系统的输出。但如果我们知道窄频带系统的输出及其两种系统的传递函数，无法得到宽频带系统的输出。这样就使得我们在记录某种信号时采用宽频带记录，然后仿真到其他各种窄频带的记录仪器上对信号进行分析。如果已知地震仪的输出和地震仪的传递函数，我们可以求出地面运动为：（10-7）下面我们就举例给出将宽频带系统的输出仿真到窄频带输出及地面运动恢复的方法。【例10-2】设计一个Butterworth模拟宽频带带通滤波器，设计指标为：通带频率：1~40Hz，低端阻带边界0.2Hz，高端阻带边界42.5Hz，通带波纹1dB，阻带衰减20dB。再设计一个窄频带带通滤波器，设计指标为通带频率：5~7Hz，低频段过渡带宽1.5Hz，高频段过渡带宽0.5Hz，通带波纹1dB，阻带衰减20dB。假设一个信号，其中f1=6Hz,f2=10Hz,f3=19Hz。信号的采样频率为50Hz。（1）模拟输入信号通过宽频带滤波器的输出信号并与原来信号进行比较。（2）运用宽带滤波器的传递函数和输出得到输入信号，并与原输入信号进行比较（恢复地面运动）（3）模拟输入信号通过窄频带滤波器的输出信号并与原来信号进行比较。（4）将在宽频带滤波器的模拟输出信号仿真到窄频带滤波器上并与（3）的结果进行比较。%Samp10_2Fs=100;dt=1/Fs;%给定模拟输入信号采样间隔f1=6;f2=10;f3=19;%模拟输入信号的频率成分N=500;%数据点数t=[0:N-1]*dt;%模拟输入信号的时间序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t);%模拟输入信号X=fft(x);%得到输入信号的Fouier变换%设计宽带Butterworth模拟滤波器wp=[140]*2*pi;ws=[0.242.5]*2*pi;%通带和阻带截止频率Rp=1;Rs=20;%通带波纹和阻带衰减[Order,Wn]=buttord(wp,ws,Rp,Rs,'s');%求解Butterworth滤波器的最小阶数w=linspace(0,Fs/2,N/2)*2*pi;%设置绘制幅频特性的频率[b1,a1]=butter(Order,Wn,'s');%设计Butterworth带通滤波器H=freqs(b1,a1,w);%计算带通滤波器的频率响应magH=abs(H);phaH=unwrap(angle(H));%求得振幅谱和相位谱figure(1)subplot(3,1,1),plot(w/(2*pi),20*log10(magH));%绘制幅频特性xlabel('频率/Hz')ylabel('振幅/dB');ylim([-1000]);xlim([0,50]);gridontitle('宽带模拟带通滤波器');%模拟输出H=freqs(b1,a1,w);%滤波器的幅频响应H1=zeros(1,N);forii=1:N/2%构建与Fourier变换得到前后数据共轭的形式H1(ii)=H(ii);H1(N-ii)=conj(H1(ii));endX=fft(x);%将原信号进行Fourier变换Y=zeros(1,N);forii=1:NY(ii)=X(ii).*H1(ii);%按（10-4）式得到的输出endy=real(ifft(Y));%将频率域的数据通过Fourier变换转换到时间域subplot(3,1,2),plot(t,x),title('输入信号')%绘制输入信号subplot(3,1,3),plot(t,y)%绘制输出信号title('宽带滤波器的模拟输出信号')xlabel('时间/s')%恢复地面运动XX=zeros(1,N);forii=1:Nif(abs(H1(ii))>1.0e-1)%为防止位于阻带内数据不稳定XX(ii)=Y(ii)./H1(ii);%恢复地面运动公式（10-7）endendxx=real(ifft(XX));%转换到时间域figure(2)plot(t,x,t,xx,':r')legend('原始信号','恢复信号',1)%给出图例title('原始信号和恢复信号')xlabel('时间/s')%窄带滤波器的设计wp=[57]*2*pi;ws=[3.57.5]*2*pi;%窄带滤波器的通带和阻带边界频率Rp=1;Rs=20;%窄带滤波器的通带波纹和阻带衰减[Order,Wn]=buttord(wp,ws,Rp,Rs,'s');%计算窄带滤波器的最小阶数和截止频率w=linspace(0,Fs/2,N/2)*2*pi;%设置绘制窄带滤波器幅频特性的频率[b2,a2]=butter(Order,Wn,'s');%设计Butterworth窄带滤波器H=freqs(b2,a2,w);%计算Butterworth窄带滤波器的频率响应magH=abs(H);phaH=unwrap(angle(H));%计算振幅谱和相位谱figure(3)subplot(3,1,1),plot(w/(2*pi),20*log10(magH));%绘制幅频特性xlabel('频率/Hz'),ylabel('振幅/dB');ylim([-1000]);xlim([0,50])title('窄带模拟滤波器');gridon%模拟输出H=freqs(b2,a2,w);H2=zeros(1,N);forii=1:N/2H2(ii)=H(ii);H2(N-ii)=conj(H2(ii));endX=fft(x);Y1=zeros(1,N);forii=1:NY1(ii)=X(ii).*H2(ii);endy1=real(ifft(Y1));subplot(3,1,2),plot(t,x),title('输入信号')%绘制输入信号subplot(3,1,3),plot(t,y1)%绘制输出信号title('窄带滤波器的模拟输出')xlabel('时间/s')%运用宽频带仪器的输出仿真得到窄带仪器上的输出figure(4)XY=zeros(1,N);forii=1:Nif(abs(H1(ii))>1.0e-2)XY(ii)=Y(ii)*H2(ii)/H1(ii);%公式（10-6）endendxy=real(ifft(XY));%转换到时间域plot(t,y1,t,xy,':r')legend('窄带输出','仿真输出',1)xlabel('时间/s')程序的输出为图10-2~5。可见运用宽频带的仪器响应和其输出信号，可以较为精确地恢复地面运动（图10-3）；将宽带滤波器的输入仿真到窄带滤波器上，得到和窄带滤波器一样的波形（图10-5），验证了上述方法的正确性。图10-2例10-2中宽频带仪器的幅频响应、输入信号和输出信号图10-3例10-2中利用宽频带仪器的输出和仪器的传递函数得到输入信号并与原输入信号进行比较图10-4例10-2中窄带仪器的幅频响应、输入信号和输出信号图10-5例10-2中运用宽带仪器的输出、宽带和窄带仪器的传递函数得到的窄带仪器的输出（仿真），并与原窄带仪器的输出进行比较下面我们用真实的地震数据和数字滤波器方法来揭示地面运动恢复和向窄带地震仪上的仿真。【例10-3】以首都圈地震记录文件200301201653.VBB的hns台站的上下向记录为例，运用椭圆滤波器来揭示地面运动恢复和仿真的概念。宽带仪器的通带边界频率为[0.00124.8]Hz，阻带边界频率为[0.0000124.9]Hz，通带波纹为1dB，阻带衰减为50dB。相当于短周期地震仪的窄带仪器的阻带边界频率为[0.014.5]Hz,通带边界频率为[0.13.8]Hz，通带波纹为1dB，阻带衰减为20dB；相当于长周期地震仪的窄带仪器用低通滤波器来表示，阻带边界频率为0.1Hz,通带边界频率为0.02Hz，通带波纹为1dB，阻带衰减为30dB。（1）采用宽频带仪器传递函数和输出信号，恢复输入信号（地面运动）；（2）将宽频带仪器的输出仿真到短周期窄带仪器上；并与窄带仪器的输出进行比较；（3）将宽频带仪器的输出仿真到长周期窄带仪器上；并与窄带仪器的输出进行比较；程序如下：%Samp10_3loadhns.dat;%假设的地震仪记录之前的地面运动Xt=hns;%把数据赋值给变量Fs=50;%设定采样率单位(HZ)dt=1/Fs;%求采样间隔单位(s)N=length(Xt);%得到序列的长度t=[0:N-1]*dt;%时间序列%设计一个椭圆滤波器，假定为宽频带地震仪ws=[0.0000124.9]*2/Fs;wp=[0.00124.8]*2/Fs;%通带和阻带边界频率（归一化频率）Rp=1;Rs=50;Nn=512;%通带波纹和阻带衰减以及绘制频率特性的数据点数[Order,Wn]=ellipord(wp,ws,Rp,Rs);%求取数字滤波器的最小阶数和归一化截止频率[b,a]=ellip(Order,Rp,Rs,Wn);%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器figure(1)[H,f]=freqz(b,a,Nn,Fs);%按传递函数系数、数据点数和采样频率求得滤波器的频率特性subplot(2,1,1),plot(f,20*log10(abs(H)))xlabel('频率/Hz');ylabel('振幅/dB');gridon;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))xlabel('频率/Hz');ylabel('相位/^o');gridon;y=filtfilt(b,a,Xt);%在宽带滤波器上的输出figure(2)subplot(2,1,1),plot(t,Xt)xlabel('时间/s'),title('输入信号');ylabel('振幅');subplot(2,1,2),plot(t,y)xlabel('时间/s'),title('输出信号');ylabel('振幅');%已知宽频带地震仪的频率特性,恢复地面运动[H,f]=freqz(b,a,N,Fs,'whole');%得到地震仪的特性Y=fft(y);XX=zeros(1,N);forii=1:N%按（10-7）式得到地面运动的频率域表示if(H(ii)>1.0e-4)XX(ii)=Y(ii)./H(ii);endenddisp=real(ifft(XX));%得到地面运动figure(3),plot(t,Xt,t,disp,'r')legend('原始信号','恢复地面运动',1)xlabel('时间/s')ylabel('振幅');%仿真到短周期地震仪上，短周期地震仪用一个窄带椭圆滤波器来表示ws=[0.014.5]*2/Fs;wp=[0.13.8]*2/Fs;%通带和阻带边界频率（归一化频率）Rp=1;Rs=20;Nn=512;%通带波纹和阻带衰减以及绘制频率特性的数据点数[order,Wn]=ellipord(wp,ws,Rp,Rs);%求取数字滤波器的最小阶数和归一化截止频率[b,a]=ellip(order,Rp,Rs,Wn);%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器figure(4)[H1,f]=freqz(b,a,Nn,Fs);%按传递函数系数、数据点数和采样频率求得滤波器的频率特性subplot(2,1,1),plot(f,20*log10(abs(H1)))xlabel('频率/Hz');ylabel('振幅/dB');gridon;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H1)))xlabel('频率/Hz');ylabel('相位/^o');gridon;figure(5)y1=filtfilt(b,a,Xt);%在窄带滤波器上的输出[H1,f]=freqz(b,a,N,Fs,'whole');%得到地震仪的特性XX1=zeros(1,N);forii=1:N%按（10-6）式得到仿真结果if(abs(H(ii))>1.0e-4)XX1(ii)=Y(ii).*H1(ii)/H(ii);endendx1=ifft(XX1);plot(t,y1,t,real(x1),'r')%绘制输入信号legend('实际输出','仿真输出',1)xlabel('时间/s');ylabel('振幅');gridon;%仿真到长周期地震仪上，长周期地震仪用一个窄带椭圆滤波器来表示ws=0.1*2/Fs;wp=0.02*2/Fs;%通带和阻带边界频率（归一化频率）Rp=1;Rs=30;Nn=512;%通带波纹和阻带衰减以及绘制频率特性的数据点数[Order,Wn]=ellipord(wp,ws,Rp,Rs);%求取数字滤波器的最小阶数和归一化截止频率[b,a]=ellip(Order,Rp,Rs,Wn);%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器figure(6)y1=filtfilt(b,a,Xt);%在窄带滤波器上的输出[H1,f]=freqz(b,a,N,Fs,'whole');%得到地震仪的特性XX1=zeros(1,N);forii=1:Nif(abs(H(ii))>1.0e-4)%为了防止H值太小将该频率的信号放大XX1(ii)=Y(ii).*H1(ii)./H(ii);%按（10-6）式得到仿真结果endendx1=ifft(XX1);plot(t,y1,t,real(x1),'r:')%绘制输入信号legend('实际输出','仿真输出',1)xlabel('时间/s');ylabel('振幅');gridon;图10-6例10-3中宽带仪器的频率特性：上图：幅频响应；下图：相频响应图10-7例10-3中宽带仪器输出（下图）和原输入信号（上图）的比较图10-8例10-3中运用宽带仪器传递函数和宽带仪器的输出信号恢复的地面运动和原输入信号的比较图10-9例10-3中短周期窄带仪器的频率特性：上图：幅频响应；下图：相频响应图10-10例10-3中运用宽带和短周期窄带仪器传递函数及宽带仪器输出仿真到窄带仪器上并与窄带仪器的实际输出进行比较图10-11例10-3中运用宽带和长周期窄带仪器传递函数及宽带仪器输出仿真到窄带仪器上并与窄带仪器的实际输出进行比较图10-6为宽带仪器的频率响应，我们设计了较宽的频带，几乎包含了输入信号的所有频率。输入信号经过宽带仪器输出后波形基本未变，但幅值发生了一些改变（图10-7）。运用宽带仪器的输出及传递函数恢复的地面运动与原输入信号相差很小（图10-8）。图10-9为相对的短周期窄带仪器的频率特性，运用宽带仪器输出信号、宽带和窄带仪器的传递函数仿真到窄带仪器上与原信号经窄带仪器的滤波结果相差很小（图10-10），由于是短周期仪器，因此仿真结果中短周期体波成分较为丰富，几乎看不到面波成分。运用宽带仪器输出信号、宽带和窄带长周期仪器的传递函数仿真到窄带仪器上与原信号经窄带仪器的滤波结果也相差不大（图10-11），由于是长周期仪器，因此仿真结果中短周期体波成分很小，长周期面波成分较为丰富。小波分析举例小波变换最早是1984年研究地震资料而使用的。由于小波分析具有能够根据分析对象自动调整有关参数的“自适应性”和能够根据观测对象自动“调焦”的特性而广泛应用于地震震相识别、地震波形反演等研究领域。一个地震时间序列f(t)是由不同强度的大小地震构成的，每个时刻t的强度可以用下式表示：(10-8)其中为奇异性指数或标度指数。奇异性是指时的数值（特别是）决定了的速度。标度性是指若变量t变成，则和是自相似的，即：（10-9）（10-9）式常称为标度律。自相似的标度律是多种尺度（即无特征尺度）现象的特征。由于中国地震资料的时间序列含有多种不同尺度，也即包含有很多层次。虽然是杂乱无章的，但小波变换好比显微镜，调节其中的放大倍数就可以清楚地看出在各个层次上的变化趋势。由于未对地震资料作谱分析，所以这是地震资料的较为真实的物理分析。10.3.1小波变换与突变一个时间序列信号的Fourier变换，是从频率轴（或波数域）来分析该信号是由哪些频率（或波数）的波组成的。地震资料组成的序列并不是平稳的，而是随时间不断变化的，此时Fourier变换就不能提取中的奇异性和突变点的信息，它只是将这些信息铺开到整个频率轴上。突变点对于很多问题是很重要的，例如，全球气候变暖从何时开始，地震活动性何时开始表现活跃等。但是小波变换Tg是将分解成具有局部特性的小波:（10-10）其中小波是将具有局部特性的小波函数通过平移和放大（放大倍数为1/a）而构成的。参数a具有时间的量纲，也称为小波尺度。设为具有低通性质的平滑函数，以它的一阶、二阶导数作为小波对f(t)作小波变换，可以证明:(10-11)(10-12)式中，表示卷积。也就是说，用、对f(t)作小波变换分别相当于f(t)被平滑后在对t求一阶或二阶导数。因此，对某一固定a值，f(t)的拐点即是的极值点，又是的过零点，由此可检测信号的急剧变化之处。本研究中所采用低通性质的平滑函数为高斯函数，根据它的一阶导数、二阶导数(墨西哥草帽小波函数)作为小波基函数进行突变点分析。各小波函数的表达式为：（10-13）之所以选择此函数是因为它对称、可微、可积，时频两域都是高斯型且呈平方型指数衰减特性，在时频两域均具有很好的局域性。该函数的MATLAB程序如下：function[x,y,z]=gdttest(a,t0,si)%wavelettransformMexscroHatt=(1:1:si);%给出时间序列x=(1-((t-t0)/a).^2).*exp(-((t-t0)/a).^2/2)/a;%二阶导数y=exp(-((t-t0)/a).^2/2)/a;%高斯函数z=(t-t0).*exp(-((t-t0)/a).^2/2)/a;%一阶导数为了形象地表明该研究所用的小波基，我们用500秒长，在250秒处，尺度为10的小波基及其1阶导数、2阶导数下面的程序来绘出：%Samp10_4si=500;t0=250;a=10;t=(1:1:si);%给出时间序列x=(1-((t-t0)/a).^2).*exp(-((t-t0)/a).^2/2)/a;%二阶导数y=exp(-((t-t0)/a).^2/2)/a;%高斯函数z=(t-t0).*exp(-((t-t0)/a).^2/2)/a;%一阶导数subplot(3,1,1),plot(t,x)title('高斯小波基')subplot(3,1,2),plot(t,y)title('小波基的一阶导数')subplot(3,1,3),plot(t,z)title('小波基的二阶导数')xlabel('时间/s')程序运行结果为图10-12，清楚地显示高斯小波基及其一阶导数、二阶导数的形状。图10-12高斯小波基（上图）、高斯小波基的一阶导数（中图）和二阶导数（下图）图10-13运用小波变换方法寻找突变点一个例子。图10-13上图为试验用原始信号，图10-13中图为运用高斯函数作为基函数的小波变换，图10-13下图为运用高斯函数的一阶导数作为基函数的小波变换。可见运用小波变换可准确地识别出突变点。因此运用高斯函数作为小波基的小波变换相当于对原来的信号做平滑处理，而运用高斯函数的一阶导数作为小波基的小波变换可精确地找到信号所对应的突变点，若脉冲向上则表示信号增大，反之表示信号减弱。其源程序如下：%Samp10_5x=zeros(1,500);%设置空矢量%给出检测数据fort=1:1:500if(t<200)x(t)=50.*exp(t/300).*sin(2*pi*0.01*t);elseif(t>=200)&(t<300)x(t)=50.*exp(t/300).*sin(2*pi*0.01*200)+30;elseif((t>=300)&(t<400))x(t)=50.*exp(t/300).*sin(2*pi*0.01*200)+10;elsex(t)=50.*exp(t/300).*sin(2*pi*0.01*200)+10+250*sin(2*pi*0.003*t);endend%结束检测数据x=x';%将检测数据转置以备后面应用x=x-mean(x);%去掉平均值a=1;%给出小波尺度为1forn=(1:1:500)t0=n*1;[g,fai,v]=gdttest(a,t0,500);%调用小波函数的计算xx(n)=sum(x.*g');%一阶导数的计算（墨西哥草帽波）yyy(n)=sum(x.*(fai)');%高斯函数的计算zzz(n)=sum(x.*v');%二阶导数的计算endsubplot(3,1,1),plot(x);%绘出原始数据波形ylabel('检测数据')subplot(3,1,2),plot(yyy);%绘出高斯函数为基函数的小波变换ylabel('高斯函数结果')subplot(3,1,3),plot(zzz);%绘出高斯函数的一阶导数为基函数的小波变换ylabel('一阶导数结果')xlabel('时间/s')图10-13利用小波变换检测突变点试验实例。上图：试验用的原数字信号。中图：以高斯函数作为基函数的小波变换。下图：以高斯函数的一阶导数作为基函数的小波变换10.3.2资料中国具有较为完整的历史地震目录。一个地区的地震能量释放显然是局部状态的一种特征。对于同样的目的，比能量释放更为常用的是应变释放。由于地震辐射能量的平方根与这次地震的应变释放成正比，因此可以用地震辐射能量的平方根表征地震的能量释放。Benioff首先使用这一概念绘制一给定地区（可以说是全球的）在给定时间内的积累应变释放曲线，并就世界浅源地震绘制了典型的应变释放曲线，又称为本尼奥夫曲线。本研究本尼奥夫应变的计算采用时振梁等的公式，表示如下（10-14）其中，为以t时刻为中心的0.1年之内的本尼奥夫应变积累，为某一地区0.1年内所有M≥5地震对应变释放的贡献；M为震级。考虑到目录的完整性，本研究对中国大陆本世纪地震目录采用M≥5的地震，地震总数2042个。对于华北地区，由于明朝初年（十四世纪后半段）政府决定由编写州志扩大到建立地方县志，各地地震灾情记录较全，故选择1500年以后M≥5地震目录，地震总数404个。由第一节的分析可知，以高斯函数为小波基对历史地震0.1年之内的本尼奥夫应变积累的小波变换对应于地震活动的平滑结果，而运用高斯函数的一阶导数为小波基的小波变换就对应于地震活动的突变点。10.3.3本尼奥夫应变资料的小波分析11900~2001年中国大陆地震资料的小波分析由于本尼奥夫应变不能明显表现为地震活动的平静期和活跃期，我们对由中国大陆M≥5地震目录得到的0.1年积累的本尼奥夫应变（图10-14a的最上一道信号），以高斯函数及其一阶导数为小波基进行小波变换。图10-14(a)得到的小波变换相当于以高斯函数作为脉冲响应的滤波器输出。随着小波尺度的增大，本尼奥夫曲线逐渐平滑，尺度为0.5年和1年的以高斯函数为小波基的小波变换可以明显分为四个地震活跃期，这与张国民所给出的年代一致（图10-14）。图10-14（b）为采用高斯函数一阶导数作为基小波对地震0.1年积累的本尼奥夫应变进行不同尺度的小波变换。纵观来看，图10-14（b）由如下特点：（1）随着尺度的增加，极值点数有多到少，并