1、 第10章數(shù)字信號(hào)處理的幾個(gè)前沿課題前面介紹了數(shù)字信號(hào)處理的基本知識(shí)，本章我們將介紹時(shí)譜分析、小波變換、地震觀測(cè)系統(tǒng)仿真與地面運(yùn)動(dòng)恢復(fù)等幾個(gè)數(shù)字信號(hào)前沿課題，以便大家在實(shí)際工作中參考。10.1時(shí)譜(倒譜)分析時(shí)譜分析(Cepstrumanalysis)是一種非線性信號(hào)處理技術(shù)，它在語(yǔ)言、圖像、和噪聲處理領(lǐng)域中都有廣泛的應(yīng)用。時(shí)譜可分為兩類(lèi)：復(fù)時(shí)譜和功率時(shí)譜。MATLAB信號(hào)處理工具箱提供復(fù)時(shí)譜分析的工具函數(shù)。復(fù)時(shí)譜(Complexcepstrum)的定義為：10-1)x(n)=1卜tnlx(ej)4d2兀-兀由上式可見(jiàn)，復(fù)時(shí)譜實(shí)際上是序列x(n)的Fourier變換取自然對(duì)數(shù)，再取Fourie

2、r逆變換，得到的復(fù)時(shí)譜仍然是一個(gè)序列。也就是說(shuō)，復(fù)時(shí)譜是x(n)從時(shí)間域至頻率域、頻率域至頻率域、頻率域至?xí)r間域的三次變換。MATLAB信號(hào)處理工具箱函數(shù)cceps用于估計(jì)一個(gè)序列x的復(fù)時(shí)譜，調(diào)用格式為:xhat=cceps(x)式中，x為輸入序列(實(shí)序列)；xhat為復(fù)時(shí)譜(復(fù)序列)。MATLAB信號(hào)處理工具箱還提供了序列實(shí)時(shí)(倒)譜的計(jì)算程序rceps，調(diào)用格式為Y=rceps(x)，其中x為實(shí)序列；y為實(shí)時(shí)譜，執(zhí)行的操作為:10-2)Cx二呵X(ej)悶由此可知，我們不能從序列x的實(shí)時(shí)譜重構(gòu)原始序列，因?yàn)閷?shí)時(shí)譜是根據(jù)序列Fourier變換的幅值計(jì)算的，丟失了相位方面的信息。但如果需要，可

3、采用最小相位模式估計(jì)原始序列。由于復(fù)時(shí)譜從復(fù)頻譜計(jì)算得到，不損失相位信息，因此復(fù)時(shí)譜是可逆的，實(shí)時(shí)譜過(guò)程是不可逆的。時(shí)譜分析技術(shù)廣泛地應(yīng)用于語(yǔ)言信號(hào)分析、同態(tài)濾波技術(shù)中。這里舉一個(gè)說(shuō)明復(fù)時(shí)譜在具有回聲信號(hào)測(cè)量中的應(yīng)用?！纠?0-1】設(shè)原信號(hào)是一個(gè)45Hz的正弦波，在傳播過(guò)程中遇到障礙產(chǎn)生回聲，回聲振幅衰減為原信號(hào)的0.5,并與原信號(hào)有0.2s的延遲。在某測(cè)點(diǎn)測(cè)到的信號(hào)是原信號(hào)和回聲信號(hào)的疊加。試使用復(fù)時(shí)譜分析該測(cè)點(diǎn)的信號(hào)。%Samp10_1t=0:0.01:1.49;%時(shí)間信號(hào)，采樣間隔為0.01ssig=sin(2*pi*45*t);%原始信號(hào)echo=0.5*zeros(1,20),sig

4、(1:length(t)-20);%在信號(hào)前面補(bǔ)加20個(gè)零，并使其振幅衰減一半作為回聲信號(hào)sigecho=sig+echo;c=cceps(sigecho);%求原始信號(hào)與回聲信號(hào)疊加的復(fù)時(shí)譜subplot(2,1,1),plot(t,sigecho);%繪制原始信號(hào)與回聲的疊加信號(hào)xlabel(時(shí)間/s);axis(O1.5Tl);gridontitle(時(shí)間域信號(hào))subplot(2,1,2),plot(t,c);%繪制其復(fù)時(shí)譜xlabel(時(shí)間/s);axis(01.5T1);gridontitle(時(shí)譜域信號(hào))程序運(yùn)行結(jié)果為圖10-1。可以看到，信號(hào)的復(fù)時(shí)譜在t=0.2秒處有一個(gè)峰值，

5、這就是回聲信號(hào)。表示運(yùn)用復(fù)時(shí)譜可以明確檢測(cè)到回聲信號(hào)。在數(shù)字地震儀監(jiān)測(cè)核爆破的工作中，如果采用兩次核爆在幾乎同一時(shí)間、同一地點(diǎn)發(fā)生，爆破產(chǎn)生的地震波經(jīng)過(guò)地球內(nèi)部傳到同一地震臺(tái)上的傳播路徑幾乎一致，并且爆破源產(chǎn)生的波形相差不大，因此在地震臺(tái)上表現(xiàn)為順序到達(dá)的兩個(gè)波形類(lèi)似的波，若后面的波形被淹沒(méi)，可以采用復(fù)時(shí)譜的方法進(jìn)行檢測(cè)，明確給出是否是兩次爆破事件，如果是兩次爆破事件，就可以明確給出兩次事件的時(shí)間間隔。時(shí)間域信號(hào)時(shí)譜域信號(hào)圖10-1正弦波與回聲信號(hào)疊加的波形和時(shí)譜形狀地震觀測(cè)系統(tǒng)的仿真和地面運(yùn)動(dòng)的恢復(fù)由前面的線性系統(tǒng)的講解可知，一個(gè)系統(tǒng)可以用它的系統(tǒng)函數(shù)或脈沖響應(yīng)來(lái)表示：y（t）=h（t）x（

6、t）（10-3）式中，x（t）為輸入信號(hào)，對(duì)于地震觀測(cè)系統(tǒng)來(lái)講為地面運(yùn)動(dòng)；y（t）為系統(tǒng)的輸出，對(duì)于地震觀測(cè)系統(tǒng)來(lái)講為地震記錄；h（t）為系統(tǒng)的脈沖響應(yīng)。在頻率域內(nèi)，根據(jù)卷積定理，該式可以表示為：Y）二H（）X（）（10-4）式中，H（&）為系統(tǒng)的傳遞函數(shù)，X6）和YC）分別為x（t）和y（t）的Fourier變換?，F(xiàn)在設(shè)想一個(gè)頻帶范圍很寬的線性系統(tǒng)，比如寬頻帶地震儀，其系統(tǒng)函數(shù)為H6）；另一個(gè)頻帶較窄的系統(tǒng)，比如短周期地震儀，其系統(tǒng)函數(shù)為H6）,對(duì)于同樣的輸入X6）有：1Y（）二H（）X（），Y（）=H（）X（）（10-5）11式中，Yi為頻帶較窄的系統(tǒng)記錄的頻譜；H16）為頻帶較窄系統(tǒng)的傳

7、遞函數(shù)。比較前面兩式可得：Y6)=1H6)X6)h()10-6)將上式變換到時(shí)間域就得到頻帶較窄系統(tǒng)的輸出y1（t）。也就是說(shuō)，如果知道了寬頻帶和窄頻帶系統(tǒng)的傳遞函數(shù)H6）和H6）,原則上可以從寬頻帶系統(tǒng)的輸出推測(cè)出頻帶較1窄系統(tǒng)的輸出。但如果我們知道窄頻帶系統(tǒng)的輸出及其兩種系統(tǒng)的傳遞函數(shù)，無(wú)法得到寬頻帶系統(tǒng)的輸出。這樣就使得我們?cè)谟涗浤撤N信號(hào)時(shí)采用寬頻帶記錄，然后仿真到其他各種窄頻帶的記錄儀器上對(duì)信號(hào)進(jìn)行分析。如果已知地震儀的輸出和地震儀的傳遞函數(shù)，我們可以求出地面運(yùn)動(dòng)為：X（）=Y（叫H（10-7）下面我們就舉例給出將寬頻帶系統(tǒng)的輸出仿真到窄頻帶輸出及地面運(yùn)動(dòng)恢復(fù)的方法?！纠?0-2】設(shè)計(jì)

8、一個(gè)Butterworth模擬寬頻帶帶通濾波器，設(shè)計(jì)指標(biāo)為：通帶頻率：140Hz，低端阻帶邊界0.2Hz，高端阻帶邊界42.5Hz，通帶波紋1dB，阻帶衰減20dB。再設(shè)計(jì)一個(gè)窄頻帶帶通濾波器，設(shè)計(jì)指標(biāo)為通帶頻率：57Hz，低頻段過(guò)渡帶寬1.5Hz，高頻段過(guò)渡帶寬0.5Hz，通帶波紋1dB，阻帶衰減20dB。假設(shè)一個(gè)信號(hào)x（t）二sin2f+0.5cos2f+0.5sin2吋t，其中f1=6Hz,f2=10Hz,f3=19Hz。信號(hào)的采123樣頻率為50Hz。（1）模擬輸入信號(hào)通過(guò)寬頻帶濾波器的輸出信號(hào)并與原來(lái)信號(hào)進(jìn)行比較。（2）運(yùn)用寬帶濾波器的傳遞函數(shù)和輸出得到輸入信號(hào)，并與原輸入信號(hào)進(jìn)行比

9、較（恢復(fù)地面運(yùn)動(dòng)）（3）模擬輸入信號(hào)通過(guò)窄頻帶濾波器的輸出信號(hào)并與原來(lái)信號(hào)進(jìn)行比較。（4）將在寬頻帶濾波器的模擬輸出信號(hào)仿真到窄頻帶濾波器上并與(3)的結(jié)果進(jìn)行比較。%Samp10_2Fs=100;dt=1/Fs;%給定模擬輸入信號(hào)采樣間隔f1=6;f2=10;f3=19;%模擬輸入信號(hào)的頻率成分N=500;%數(shù)據(jù)點(diǎn)數(shù)t=0:N-1*dt;%模擬輸入信號(hào)的時(shí)間序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t);%模擬輸入信號(hào)X=fft(x);%得到輸入信號(hào)的Fouier變換%設(shè)計(jì)寬帶Butterworth模擬濾波器wp=140*2*

10、pi;ws=0.242.5*2*pi;%通帶和阻帶截止頻率Rp=1;Rs=20;%通帶波紋和阻帶衰減Order,Wn=buttord(wp,ws,Rp,Rs,s);%求解Butterworth濾波器的最小階數(shù)w=linspace(0,Fs/2,N/2)*2*pi;%設(shè)置繪制幅頻特性的頻率b1,a1=butter(Order,Wn,s);%設(shè)計(jì)Butterworth帶通濾波器H=freqs(b1,a1,w);%計(jì)算帶通濾波器的頻率響應(yīng)magH=abs(H);phaH=unwrap(angle(H);%求得振幅譜和相位譜figure(1)subplot(3,1,1),plot(w/(2*pi),2

11、0*log10(magH);%繪制幅頻特性xlabel(頻率/Hz)ylabel(振幅/dB);ylim(-1000);xlim(0,50);gridontitle(寬帶模擬帶通濾波器);%模擬輸出H=freqs(b1,a1,w);%濾波器的幅頻響應(yīng)H1=zeros(1,N);forii=1:N/2%構(gòu)建與Fourier變換得到前后數(shù)據(jù)共軛的形式H1(ii)=H(ii);H1(N-ii)=conj(H1(ii);endX=fft(x);%將原信號(hào)進(jìn)行Fourier變換Y=zeros(1,N);forii=1:NY(ii)=X(ii).*H1(ii);%按(10-4)式得到的輸出endy=rea

12、l(ifft(Y);%將頻率域的數(shù)據(jù)通過(guò)Fourier變換轉(zhuǎn)換到時(shí)間域subplot(3,l,2),plot(t,x),title(輸入信號(hào))繪制輸入信號(hào)subplot(3,1,3),plot(t,y)%繪制輸出信號(hào)title(寬帶濾波器的模擬輸出信號(hào))xlabel(時(shí)間/s)%恢復(fù)地面運(yùn)動(dòng)XX=zeros(1,N);forii=1:Nif(abs(H1(ii)1.0e-1)%為防止位于阻帶內(nèi)數(shù)據(jù)不穩(wěn)定XX(ii)=Y(ii)./H1(ii);%恢復(fù)地面運(yùn)動(dòng)公式(10-7)endendxx=real(ifft(XX);%轉(zhuǎn)換到時(shí)間域figure(2)plot(t,x,t,xx,:r)legen

13、d(原始信號(hào)，恢復(fù)信號(hào)，1)%給出圖例title(原始信號(hào)和恢復(fù)信號(hào))xlabel(時(shí)間/s)%窄帶濾波器的設(shè)計(jì)wp=57*2*pi;ws=3.57.5*2*pi;%窄帶濾波器的通帶和阻帶邊界頻率Rp=1;Rs=20;%窄帶濾波器的通帶波紋和阻帶衰減Order，Wn=buttord(wp，ws，Rp，Rs，s);%計(jì)算窄帶濾波器的最小階數(shù)和截止頻率w=linspace(0，F(xiàn)s/2，N/2)*2*pi;%設(shè)置繪制窄帶濾波器幅頻特性的頻率b2，a2=butter(Order，Wn，s);%設(shè)計(jì)Butterworth窄帶濾波器H二freqs(b2,a2,w);%計(jì)算Butterworth窄帶濾波器

14、的頻率響應(yīng)magH=abs(H);phaH=unwrap(angle(H);%計(jì)算振幅譜和相位譜figure(3)subplot(3,1,1),plot(w/(2*pi),20*log10(magH);%繪制幅頻特性xlabel(頻率/Hz),ylabel(振幅/dB);ylim(T000);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)

15、=X(ii).*H2(ii);endy1=real(ifft(Y1);subplot(3,1,2),plot(t,x),title(輸入信號(hào))繪制輸入信號(hào)subplot(3,1,3),plot(t,y1)%繪制輸出信號(hào)title(窄帶濾波器的模擬輸出)xlabel(時(shí)間/s)%運(yùn)用寬頻帶儀器的輸出仿真得到窄帶儀器上的輸出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);%轉(zhuǎn)換到時(shí)間域plot(t,y1,t,xy,:r)leg

16、end(窄帶輸出，仿真輸出，1)xlabel(時(shí)間/s)程序的輸出為圖10-25?？梢?jiàn)運(yùn)用寬頻帶的儀器響應(yīng)和其輸出信號(hào)，可以較為精確地恢復(fù)地面運(yùn)動(dòng)(圖10-3)；將寬帶濾波器的輸入仿真到窄帶濾波器上，得到和窄帶濾波器一樣的波形(圖10-5)，驗(yàn)證了上述方法的正確性。寬帶濾波器的模擬輸出信號(hào)00.511.522.533.544.55時(shí)間/sO圖10-2例10-2中寬頻帶儀器的幅頻響應(yīng)、輸入信號(hào)和輸出信號(hào)2 原始信號(hào)和恢復(fù)信號(hào)1.510.50-0.5-1-1.5-200.511.522.533.544.55時(shí)間/s圖10-3例10-2中利用寬頻帶儀器的輸出和儀器的傳遞函數(shù)得到輸入信號(hào)并與原輸入信號(hào)

17、進(jìn)行比較 窄帶模擬濾波器輸入信號(hào)0.50IJIIIJ0.51.522.533.544.55時(shí)間/s圖10-4例10-2中窄帶儀器的幅頻響應(yīng)、輸入信號(hào)和輸出信號(hào)原窄帶儀器的輸出進(jìn)行比較下面我們用真實(shí)的地震數(shù)據(jù)和數(shù)字濾波器方法來(lái)揭示地面運(yùn)動(dòng)恢復(fù)和向窄帶地震儀上的仿真。【例10-3】以首都圈地震記錄文件200301201653.VBB的hns臺(tái)站的上下向記錄為例，運(yùn)用橢圓濾波器來(lái)揭示地面運(yùn)動(dòng)恢復(fù)和仿真的概念。寬帶儀器的通帶邊界頻率為0.00124.8Hz,阻帶邊界頻率為0.0000124.9Hz,通帶波紋為1dB,阻帶衰減為50dB。相當(dāng)于短周期地震儀的窄帶儀器的阻帶邊界頻率為0.014.5Hz,通

18、帶邊界頻率為0.13.8Hz,通帶波紋為1dB，阻帶衰減為20dB；相當(dāng)于長(zhǎng)周期地震儀的窄帶儀器用低通濾波器來(lái)表示，阻帶邊界頻率為0.1Hz,通帶邊界頻率為0.02Hz，通帶波紋為1dB，阻帶衰減為30dB。（1）采用寬頻帶儀器傳遞函數(shù)和輸出信號(hào),恢復(fù)輸入信號(hào)（地面運(yùn)動(dòng)）；（2）將寬頻帶儀器的輸出仿真到短周期窄帶儀器上；并與窄帶儀器的輸出進(jìn)行比較；（3）將寬頻帶儀器的輸出仿真到長(zhǎng)周期窄帶儀器上；并與窄帶儀器的輸出進(jìn)行比較；%Samp10_3loadhns.datXt=hns;Fs=50;程序如下：%假設(shè)的地震儀記錄之前的地面運(yùn)動(dòng)%把數(shù)據(jù)賦值給變量%設(shè)定采樣率單位（HZ）dt=1/Fs;%求采樣

19、間隔單位（s）N=length（Xt）;%得到序列的長(zhǎng)度t=0:N-1*dt;%時(shí)間序列%設(shè)計(jì)一個(gè)橢圓濾波器，假定為寬頻帶地震儀ws=0.0000124.9*2/Fs;wp=0.00124.8*2/Fs;%通帶和阻帶邊界頻率(歸一化頻率)Rp=1;Rs=50;Nn=512;%通帶波紋和阻帶衰減以及繪制頻率特性的數(shù)據(jù)點(diǎn)數(shù)Order,Wn=ellipord(wp,ws,Rp,Rs);%求取數(shù)字濾波器的最小階數(shù)和歸一化截止頻率b,a=ellip(Order,Rp,Rs,Wn);%按最小階數(shù)、截止頻率、通帶波紋和阻帶衰減設(shè)計(jì)濾波器figure(1)H,f=freqz(b,a,Nn,Fs);%按傳遞函數(shù)

20、系數(shù)、數(shù)據(jù)點(diǎn)數(shù)和采樣頻率求得濾波器的頻率特性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(時(shí)間/s),title(輸入信號(hào));ylabel(振幅);subplot(2,1,2),plot(t,y)xlab

21、el(時(shí)間/s),title(輸出信號(hào));ylabel(振幅);%已知寬頻帶地震儀的頻率特性,恢復(fù)地面運(yùn)動(dòng)H,f=freqz(b,a,N,Fs,whole);%得到地震儀的特性Y=fft(y);XX=zeros(1,N);forii=1:N%按(10-7)式得到地面運(yùn)動(dòng)的頻率域表示if(H(ii)1.0e-4)XX(ii)=Y(ii)./H(ii);endenddisp=real(ifft(XX);%得到地面運(yùn)動(dòng)figure(3),plot(t,Xt,t,disp,r)legend(原始信號(hào),恢復(fù)地面運(yùn)動(dòng),1)xlabel(時(shí)間/s)ylabel(振幅);%仿真到短周期地震儀上，短周期地震儀用

22、一個(gè)窄帶橢圓濾波器來(lái)表示ws=0.014.5*2/Fs;wp=0.13.8*2/Fs;%通帶和阻帶邊界頻率(歸一化頻率)Rp=1;Rs=20;Nn=512;%通帶波紋和阻帶衰減以及繪制頻率特性的數(shù)據(jù)點(diǎn)數(shù)order,Wn=ellipord(wp,ws,Rp,Rs);%求取數(shù)字濾波器的最小階數(shù)和歸一化截止頻率b,a=ellip(order,Rp,Rs,Wn);%按最小階數(shù)、截止頻率、通帶波紋和阻帶衰減設(shè)計(jì)濾波器figure(4)H1,f=freqz(b,a,Nn,Fs);%按傳遞函數(shù)系數(shù)、數(shù)據(jù)點(diǎn)數(shù)和采樣頻率求得濾波器的頻率特性subplot(2,1,1),plot(f,20*log10(abs(H

23、1)xlabel(頻率/Hz);ylabel(振幅/dB);gridon;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H1)xlabel(頻率/Hz);ylabel(相位,/Ao);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)式得到仿真結(jié)果if(abs(H(ii)1.0e-4)XX1(ii)=Y(ii).*H1(ii)/H(ii);endendx1=ifft(XX1);p

24、lot(t,y1,t,real(x1),r)%繪制輸入信號(hào)legend(實(shí)際輸出，仿真輸出，1)xlabel(時(shí)間/s);ylabel(振幅)；gridon;%仿真到長(zhǎng)周期地震儀上，長(zhǎng)周期地震儀用一個(gè)窄帶橢圓濾波器來(lái)表示ws=0.1*2/Fs;wp=0.02*2/Fs;%通帶和阻帶邊界頻率(歸一化頻率)Rp=1;Rs=30;Nn=512;%通帶波紋和阻帶衰減以及繪制頻率特性的數(shù)據(jù)點(diǎn)數(shù)Order，Wn=ellipord(wp，ws，Rp，Rs);%求取數(shù)字濾波器的最小階數(shù)和歸一化截止頻率b，a=ellip(Order，Rp，Rs，Wn);%按最小階數(shù)、截止頻率、通帶波紋和阻帶衰減設(shè)計(jì)濾波器fig

25、ure(6)y1=filtfilt(b，a，Xt);%在窄帶濾波器上的輸出H1，f=freqz(b，a，N，F(xiàn)s，whole);%得到地震儀的特性XX1=zeros(1，N);forii=1:Nif(abs(H(ii)1.0e-4)%為了防止H值太小將該頻率的信號(hào)放大XX1(ii)=Y(ii).*H1(ii)./H(ii);%按(10-6)式得到仿真結(jié)果endendx1=ifft(XX1);plot(t，y1，t，real(x1)，r:)%繪制輸入信號(hào)legend(實(shí)際輸出，仿真輸出，1)xlabel(時(shí)間/s);ylabel(振幅)；gridon;圖10-6例10-3中寬帶儀器的頻率特性：上

26、圖：幅頻響應(yīng)；下圖：相頻響應(yīng)642o-2-45叩I刪0幅振幅振60-50 x104輸入信號(hào)5002000250010001500時(shí)間/s5002000250010001500時(shí)間/sx104輸出信號(hào)圖10-7例10-3中寬帶儀器輸出（下圖）和原輸入信號(hào)（上圖）的比較圖10-9例10-3中短周期窄帶儀器的頻率特性：上圖：幅頻響應(yīng)；下圖：相頻響應(yīng)x104時(shí)間/s圖10-10例10-3中運(yùn)用寬帶和短周期窄帶儀器傳遞函數(shù)及寬帶儀器輸出仿真到窄帶儀器上并與窄帶儀器的實(shí)際輸出進(jìn)行比較x104時(shí)間/s圖10-11例10-3中運(yùn)用寬帶和長(zhǎng)周期窄帶儀器傳遞函數(shù)及寬帶儀器輸出仿真到窄帶儀器上并與窄帶儀器的 實(shí)際輸

27、出進(jìn)行比較圖10-6為寬帶儀器的頻率響應(yīng)，我們?cè)O(shè)計(jì)了較寬的頻帶，幾乎包含了輸入信號(hào)的所有頻率。輸入信號(hào)經(jīng)過(guò)寬帶儀器輸出后波形基本未變，但幅值發(fā)生了一些改變(圖10-7)。運(yùn)用寬帶儀器的輸出及傳遞函數(shù)恢復(fù)的地面運(yùn)動(dòng)與原輸入信號(hào)相差很小(圖10-8)。圖10-9為相對(duì)的短周期窄帶儀器的頻率特性，運(yùn)用寬帶儀器輸出信號(hào)、寬帶和窄帶儀器的傳遞函數(shù)仿真到窄帶儀器上與原信號(hào)經(jīng)窄帶儀器的濾波結(jié)果相差很小(圖10-10)，由于是短周期儀器，因此仿真結(jié)果中短周期體波成分較為豐富，幾乎看不到面波成分。運(yùn)用寬帶儀器輸出信號(hào)、寬帶和窄帶長(zhǎng)周期儀器的傳遞函數(shù)仿真到窄帶儀器上與原信號(hào)經(jīng)窄帶儀器的濾波結(jié)果也相差不大(圖10-

28、11)，由于是長(zhǎng)周期儀器，因此仿真結(jié)果中短周期體波成分很小，長(zhǎng)周期面波成分較為豐富。10.3小波分析舉例小波變換最早是1984年研究地震資料而使用的。由于小波分析具有能夠根據(jù)分析對(duì)象自動(dòng)調(diào)整有關(guān)參數(shù)的“自適應(yīng)性”和能夠根據(jù)觀測(cè)對(duì)象自動(dòng)“調(diào)焦”的特性而廣泛應(yīng)用于地震震相識(shí)別、地震波形反演等研究領(lǐng)域。一個(gè)地震時(shí)間序列f(t)是由不同強(qiáng)度的大小地震構(gòu)成的，每個(gè)時(shí)刻t的強(qiáng)度可以用下式表示：f(t)tQ(10-8)其中Q為奇異性指數(shù)或標(biāo)度指數(shù)。奇異性是指tT0時(shí)Q的數(shù)值(特別是Q0)決定了f(t)T0的速度。標(biāo)度性是指若變量t變成九t，則f(t)和f(九t)是自相似的，即：f(kt)二(XtL二Aqta

29、二辰f(t)(10-9)(10-9)式常稱(chēng)為標(biāo)度律。自相似的標(biāo)度律是多種尺度(即無(wú)特征尺度)現(xiàn)象的特征。由于中國(guó)地震資料的時(shí)間序列含有多種不同尺度，也即包含有很多層次。雖然f(t)是雜亂無(wú)章的，但小波變換好比顯微鏡，調(diào)節(jié)其中的放大倍數(shù)就可以清楚地看出f(t)在各個(gè)層次上的變化趨勢(shì)。由于未對(duì)地震資料作譜分析，所以這是地震資料的較為真實(shí)的物理分析。10.3.1小波變換與突變一個(gè)時(shí)間序列信號(hào)f(t)的Fourier變換，是從頻率軸(或波數(shù)域)來(lái)分析該信號(hào)是由哪些頻率(或波數(shù))的波組成的。地震資料f(t)組成的序列并不是平穩(wěn)的，而是隨時(shí)間不斷變化的，此時(shí)Fourier變換就不能提取f(t)中的奇異性和

30、突變點(diǎn)的信息，它只是將這些信息鋪開(kāi)到整個(gè)頻率軸上。突變點(diǎn)對(duì)于很多問(wèn)題是很重要的，例如，全球氣候變暖從何時(shí)開(kāi)始， 地震活動(dòng)性何時(shí)開(kāi)始表現(xiàn)活躍等。但是小波變換Tg是將f(t)分解成具有局部特性的小波gQ：T(a,t)=g0If(tt)-f(t)gdtag.a丿-0-0)-(其中小波一gatt0aI是將具有局部特性的小波函數(shù)g通過(guò)平移和放大(放大倍數(shù)為1/a)而構(gòu)成的。參數(shù)a具有時(shí)間的量綱，也稱(chēng)為小波尺度。設(shè)就)為具有低通性質(zhì)的平滑函數(shù)，以它的一階、二階導(dǎo)數(shù)作為小波對(duì)f(t)作小波變換,可以證明:d(10-11)Wf(a,T)二f(t)-(at)二a-f(t)6(t)TOC o 1-5 h z91d

31、tWf(a,T)=f(t)申2(at)=af(t)00(t)(10-12) HYPERLINK l bookmark64 申2dt2式中，表示卷積。也就是說(shuō)，用1(t)、cp2(t)對(duì)f(t)作小波變換分別相當(dāng)于f(t)被0C)平滑后在對(duì)t求一階或二階導(dǎo)數(shù)。因此，對(duì)某一固定a值,f(t)00C)的拐點(diǎn)即是Wf(a,T)ap1的極值點(diǎn)，又是Wf(a,T)的過(guò)零點(diǎn)，由此可檢測(cè)信號(hào)的急劇變化之處。p2本研究中所采用低通性質(zhì)的平滑函數(shù)為高斯函數(shù)，根據(jù)它的一階導(dǎo)數(shù)、二階導(dǎo)數(shù)(墨西哥草帽小波函數(shù))作為小波基函數(shù)進(jìn)行突變點(diǎn)分析。各小波函數(shù)的表達(dá)式為：t王10-13)g(t,a)=e2a2a之所以選擇此函數(shù)是

32、因?yàn)樗鼘?duì)稱(chēng)、可微、可積，時(shí)頻兩域都是高斯型且呈平方型指數(shù)衰減特性，在時(shí)頻兩域均具有很好的局域性。該函數(shù)的MATLAB程序如下：functionx,y,z=gdttest(a,t0,si)%wavelettransformMexscroHatt=(1:1:si);%給出時(shí)間序列x=(l-(t-tO)/a).2).*exp(-(t-t0)/a).2/2)/a;%二階導(dǎo)數(shù)y=exp(-(t-tO)/a).2/2)/a;%高斯函數(shù)z=(t-tO).*exp(-(t-tO)/a).2/2)/a;%一階導(dǎo)數(shù)為了形象地表明該研究所用的小波基，我們用500秒長(zhǎng)，在250秒處，尺度為10的小波基及其1階導(dǎo)數(shù)、2

33、階導(dǎo)數(shù)下面的程序來(lái)繪出：%Samp10_4si=500;t0=250;a=10;t=(1:1:si);%給出時(shí)間序列x=(l-(t-tO)/a).2).*exp(-(t-t0)/a).2/2)/a;%二階導(dǎo)數(shù)y=exp(-(t-tO)/a).2/2)/a;%高斯函數(shù)z=(t-tO).*exp(-(t-tO)/a).2/2)/a;%一階導(dǎo)數(shù)subplot(3,1,1),plot(t,x)title(高斯小波基)subplot(3,1,2),plot(t,y)title(小波基的一階導(dǎo)數(shù))subplot(3,1,3),plot(t,z)title(小波基的二階導(dǎo)數(shù))xlabel(時(shí)間/s)程序運(yùn)行

34、結(jié)果為圖10-12，清楚地顯示高斯小波基及其一階導(dǎo)數(shù)、二階導(dǎo)數(shù)的形狀。-0.05咼斯小波基0.150.10.050時(shí)間/s圖10-12咼斯小波基（上圖）、咼斯小波基的一階導(dǎo)數(shù)（中圖）和二階導(dǎo)數(shù)（下圖）圖10-13運(yùn)用小波變換方法尋找突變點(diǎn)一個(gè)例子。圖10-13上圖為試驗(yàn)用原始信號(hào)，圖10-13中圖為運(yùn)用咼斯函數(shù)作為基函數(shù)的小波變換，圖10-13下圖為運(yùn)用咼斯函數(shù)的一階導(dǎo)數(shù)作為基函數(shù)的小波變換?？梢?jiàn)運(yùn)用小波變換可準(zhǔn)確地識(shí)別出突變點(diǎn)。因此運(yùn)用咼斯函數(shù)作為小波基的小波變換相當(dāng)于對(duì)原來(lái)的信號(hào)做平滑處理，而運(yùn)用咼斯函數(shù)的一階導(dǎo)數(shù)作為小波基的小波變換可精確地找到信號(hào)所對(duì)應(yīng)的突變點(diǎn)，若脈沖向上則表示信號(hào)增大

35、，反之表示信號(hào)減弱。其源程序如下：%Samp10_5x=zeros（1,500）;%設(shè)置空矢量%給出檢測(cè)數(shù)據(jù)fort=1:1:500if(t=200)&(t=300)&(t400)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%結(jié)束檢測(cè)數(shù)據(jù)x=x;%將檢測(cè)數(shù)據(jù)轉(zhuǎn)置以備后面應(yīng)用x=x-mean(x);%去掉平均值a=1;%給出小波尺度為1forn=(1:1:500)t0=n*1;g,fai,v=gdttest

36、(a,t0,500);%調(diào)用小波函數(shù)的計(jì)算xx(n)=sum(x.*g);%一階導(dǎo)數(shù)的計(jì)算(墨西哥草帽波)yyy(n)=sum(x.*(fai);%高斯函數(shù)的計(jì)算zzz(n)=sum(x.*v);%二階導(dǎo)數(shù)的計(jì)算endsubplot(3,1,1),plot(x);%繪出原始數(shù)據(jù)波形ylabel(檢測(cè)數(shù)據(jù))subplot(3,1,2),plot(yyy);%繪出高斯函數(shù)為基函數(shù)的小波變換ylabel(高斯函數(shù)結(jié)果)subplot(3,1,3),plot(zzz);%繪出高斯函數(shù)的一階導(dǎo)數(shù)為基函數(shù)的小波變換ylabel(階導(dǎo)數(shù)結(jié)果)xlabel(時(shí)間/s)據(jù)數(shù)測(cè)檢果結(jié)數(shù)函斯高oooOooO321果

37、結(jié)數(shù)導(dǎo)階一050100150200250300350400450500時(shí)間/s圖10-13利用小波變換檢測(cè)突變點(diǎn)試驗(yàn)實(shí)例。上圖：試驗(yàn)用的原數(shù)字信號(hào)。中圖：以高斯函數(shù)作為基函數(shù)的小波變換。下圖：以高斯函數(shù)的一階導(dǎo)數(shù)作為基函數(shù)的小波變換資料中國(guó)具有較為完整的歷史地震目錄。一個(gè)地區(qū)的地震能量釋放顯然是局部狀態(tài)的一種特征。對(duì)于同樣的目的，比能量釋放更為常用的是應(yīng)變釋放。由于地震輻射能量的平方根與這次地震的應(yīng)變釋放成正比，因此可以用地震輻射能量的平方根表征地震的能量釋放。Benioff首先使用這一概念繪制一給定地區(qū)（可以說(shuō)是全球的）在給定時(shí)間內(nèi)的積累應(yīng)變釋放曲線，并就世界淺源地震繪制了典型的應(yīng)變釋放曲線

38、，又稱(chēng)為本尼奧夫曲線。本研究本尼奧夫應(yīng)變的計(jì)算采用時(shí)振梁等的公式，表示如下（t）=Z10（ZM）2（10-14）b其中，e（t）為以t時(shí)刻為中心的0.1年之內(nèi)的本尼奧夫應(yīng)變積累，為某一地區(qū)0.1年內(nèi)b所有M25地震對(duì)應(yīng)變釋放的貢獻(xiàn)；M為震級(jí)?？紤]到目錄的完整性，本研究對(duì)中國(guó)大陸本世紀(jì)地震目錄采用M25的地震，地震總數(shù)2042個(gè)。對(duì)于華北地區(qū)，由于明朝初年（十四世紀(jì)后半段）政府決定由編寫(xiě)州志擴(kuò)大到建立地方縣志，各地地震災(zāi)情記錄較全，故選擇1500年以后M25地震目錄，地震總數(shù)404個(gè)。由第一節(jié)的分析可知，以高斯函數(shù)為小波基對(duì)歷史地震0.1年之內(nèi)的本尼奧夫應(yīng)變積累的小波變換對(duì)應(yīng)于地震活動(dòng)的平滑結(jié)果

39、，而運(yùn)用高斯函數(shù)的一階導(dǎo)數(shù)為小波基的小波變換就對(duì)應(yīng)于地震活動(dòng)的突變點(diǎn)。本尼奧夫應(yīng)變資料的小波分析119002001年中國(guó)大陸地震資料的小波分析由于本尼奧夫應(yīng)變不能明顯表現(xiàn)為地震活動(dòng)的平靜期和活躍期，我們對(duì)由中國(guó)大陸M三5地震目錄得到的0.1年積累的本尼奧夫應(yīng)變（圖10-14a的最上一道信號(hào)），以高斯函數(shù)及其一階導(dǎo)數(shù)為小波基進(jìn)行小波變換。圖10-14（a）得到的小波變換相當(dāng)于以高斯函數(shù)作為脈沖響應(yīng)的濾波器輸出。隨著小波尺度的增大，本尼奧夫曲線逐漸平滑，尺度為0.5年和1年的以高斯函數(shù)為小波基的小波變換可以明顯分為四個(gè)地震活躍期，這與張國(guó)民所給出的年代一致（圖10-14）。圖10-14（b）為采用

40、高斯函數(shù)一階導(dǎo)數(shù)作為基小波對(duì)地震0.1年積累的本尼奧夫應(yīng)變進(jìn)行不同尺度的小波變換?？v觀來(lái)看，圖10-14（b）由如下特點(diǎn)：（1）隨著尺度的增加，極值點(diǎn)數(shù)有多到少，并且相鄰的極值點(diǎn)逐漸融合成較大的極值點(diǎn)（。2）各極值點(diǎn)對(duì)應(yīng)于運(yùn)用高斯函數(shù)的相同尺度進(jìn)行的小波變換（濾波后的本尼奧夫曲線，圖10-14a）上的局部最大突變點(diǎn)（即原信號(hào)的拐點(diǎn)），極值點(diǎn)的幅值反映了對(duì)應(yīng)點(diǎn)局部梯度的大小，極值點(diǎn)的正負(fù)反映了濾波后本尼奧夫曲線的邊沿方向特性，正極大值點(diǎn)對(duì)應(yīng)于的上升沿，負(fù)極大值點(diǎn)對(duì)應(yīng)于下降沿。（3）極值點(diǎn)隨尺度變化的幅值大小反映了曲線邊沿變化特性，幅值較大反映了濾波后本尼奧夫曲線較劇烈的變化趨勢(shì)（較陡的邊沿），較小的幅值對(duì)應(yīng)于濾波后本尼奧夫曲線較平緩的變化趨勢(shì)（較平的邊沿）。由此可得到地震活躍期的開(kāi)始年代（地震突變點(diǎn)）（4）當(dāng)尺度增大時(shí)極值點(diǎn)消失或成對(duì)結(jié)合成新的極值點(diǎn)，反映了信號(hào)中的隱含的某種“嵌套結(jié)構(gòu)”嵌套范圍大的曲線反映了原始信號(hào)中大范圍能表現(xiàn)的特征。根據(jù)上面的分析可以看出，1900年以來(lái)中國(guó)大陸地震活動(dòng)的第一個(gè)活躍期，對(duì)應(yīng)于1902年新疆阿圖什地震和1906年新疆瑪納斯地