信號與系統(tǒng)課程設計報告-管道泄漏檢測定位應用中的信號分析與處理

1、<p>　　《信號與系統(tǒng)》課程設計報告</p><p>　　管道泄漏檢測定位應用中的信號分析與處理</p><p>　　學院： </p><p>　　班級： </p><p>　　姓名： </p><p>　　學號：

2、 </p><p>　　指導老師： </p><p><b>　　本小組其他成員： </b></p><p><b>　　摘要</b></p><p>　　管道泄漏時，泄漏點將會產生振動，出現(xiàn)聲音信號。</p><p>　　通過對兩個采樣器采集

3、到的振動信號進行時域和頻域的分析，使用巴特沃茲高通濾波器，濾掉低頻（50Hz）噪聲信號。采用互相關函數找到兩個采樣器接收信號的時間差，根據聲音在管道中的速度，求出管道泄漏點的位置。利用MATLAB實現(xiàn)以上功能。</p><p>　　關鍵字：時頻域分析、巴特沃茲高通濾波器、互相關函數</p><p><b>　　二、設計原理及要求</b></p><

4、p>　　【管道泄漏檢測原理】</p><p>　　已知：某學校有段自來水管道泄漏，該管道為鑄鐵管道，管徑為250mm。在管道兩端接頭處分別放置一個傳感器及采集模塊進行采集管道振動數據進行分析處理，得到漏點位置。已知兩個采集器的采樣率均為10KHz，測試點之間的管道總長度L=100m，水聲在管徑為250mm鑄鐵管道中的傳播速度v=1.250m/ms。試圖采用互相關方法對泄漏信號進行分析得到漏點位置，檢測與定位

5、原理如圖1：</p><p>　　圖1 管道泄漏檢測與定位原理</p><p><b>　　注：</b></p><p>　　L1, L2 = 漏點位置離傳感器1、2的距離; </p><p>　　L = 兩傳感器間的管道長度; </p><p>　　v = 泄漏噪聲在管道內的傳播速度; &l

6、t;/p><p>　　D = 泄漏聲信號到達兩傳感器的時間差. </p><p><b>　　泄漏點定位公式：</b></p><p><b>　　【設計要求】</b></p><p>　　對兩端采集的數據分別進行時頻域分析；</p><p>　　對兩端采集的數據進行互相關運算，

7、計算兩路信號的相對延遲時間，并根據已給的兩傳感器間隔距離及水聲沿管道傳播速度計算泄漏點位置。</p><p><b>　　【設計步驟及提示】</b></p><p>　　通過MATLAB的函數load ()可以讀入兩個.txt格式的數據文件，并將該文件保存到指定的數組中。例如下面的語句（更詳細的命令介紹可以自己查閱MATLAB的幫助）中，將.txt讀入后分別存放到矩陣

8、y1和y2中。</p><p>　　y1 = load(' D:\MATLAB\R2006a\work\S&Sdesign\2005年12月2日15點39分set2daq1\data1.txt'); </p><p>　　y2 = load(' D:\MATLAB\R2006a\work\S&Sdesign\2005年12月2日15點39分set2

9、daq1\data1.txt'); </p><p>　　對泄漏信號進行時頻域分析。要求編寫一個簡單的函數文件SignalTF.m，它具有如下性質：a、函數實現(xiàn)對信號進行快速傅里葉變換(FFT)；b、繪制時域波形與頻譜圖。要求圖形中坐標標注齊全：1）時域波形坐標標注為Time（s）—Amplitude（v），標題為“waveform”2）頻譜圖坐標標注為Frequency（Hz）—Magnitude（v

10、/Hz）（即，橫坐標標注為Frequency（Hz），縱坐標標注為Magnitude），標題為“Frequency spectrum”。找到泄漏信號的主要頻譜成分所在的頻帶，驗證為何檢漏系統(tǒng)可以對泄漏信號采用10KHz的采樣速率，并指出最低的采樣頻率是多少。</p><p>　　利用MATLAB函數庫中函數設計butterworth高通濾波器（見Matlab 的butter, filter函數），對采集的兩路信號

11、分別進行高通濾波，抑制工頻噪聲（典型值50Hz）及其他常見低頻干擾噪聲。要求編寫PreProcess.m函數，a) 通過畫圖對濾波前后信號進行對比，觀察差異，并用subplot命令分別將濾波前后的信號畫在一個FIGURE圖中（兩個信號均需給出其濾波前后圖形）；b) 將濾波后的兩組數據分別寫至data1p.txt和data2p.txt中。</p><p>　　編寫TDE.m函數，用互相關方法檢測兩路信號的相對時間延

12、遲，并以此計算漏點位置。（提示：用Matlab中xcorr函數實現(xiàn)兩個信號互相關運算，根據峰值位置確定兩路信號的相對時間間隔即數據點，并根據采樣頻率計算延遲時間，以此計算漏點位置。）要求：a) 畫出兩路信號的互相關曲線圖，找出最大點或峰值點及其橫坐標，并標注齊全，峰值位置需對應正確的延遲時間（提示：用Matlab中find函數查找最大值點及其位置）。</p><p><b>　　信號說明：</b&

13、gt;</p><p>　　data1.txt：傳感器1采集的泄漏信號；data2.txt：傳感器2采集的泄漏信號。</p><p>　　【MATLAB參考命令】</p><p>　　繪圖命令：plot、subplot等</p><p>　　求頻譜、頻率響應：fft、freqz等</p><p>　　卷積、濾波：con

14、v、deconv、filter等</p><p>　　文件讀寫：load、fopen, fprintf, fclose, xlsread, xlswrite等</p><p><b>　　三、正文</b></p><p>　　1.對泄露信號進行時頻域分析</p><p>　　在MATLAB中，通過stem函數，可以畫出離

15、散時間信號在時域上的圖像。通過fft函數實現(xiàn)信號的快速傅里葉變換，使得信號從時域變換到頻域上。使用fftshift函數實現(xiàn)頻率搬移，使頻譜圖關于原點左右對稱。這道題目的難點是確定頻譜圖的橫坐標（即頻率）的范圍。由采樣定理，信號的最高頻率不得超過采樣頻率的一半，否則采樣后的信號將會出現(xiàn)失真。由此可以確定頻譜圖的頻率范圍為-fs/2 到 fs/2（fs為采樣頻率）。用函數linspace(-fs/2,fs/2,num)實現(xiàn)（num為采樣信號

16、的個數） </p><p>　　得到兩個信號的時域圖像如下：</p><p>　　得到信號的頻譜圖分別如下：</p><p>　　從頻譜圖可以看出，兩個信號的最高頻率大約為2200Hz，最低采樣頻率Fsmin=2*fm=4400Hz</p><p><b>　　程序代碼如下：</b></p><p&

17、gt;　　function Y=SignalTF(X,N) %N為采樣周期</p><p>　　n=N:N:length(X)*N; %得到離散時間序列</p><p><b>　　fs=1/N;</b></p><p>　　num=length(X);</p><p><b>　　Y=fft(X);&

18、lt;/b></p><p>　　Y1=fftshift(Y);</p><p><b>　　figure,</b></p><p>　　stem(n,X,'filled') %畫出時域圖像</p><p>　　title('waveform')</p><

19、;p>　　xlabel('Time(s)')</p><p>　　ylabel('Amlitude(v)')</p><p>　　mag=abs(Y1); %幅度</p><p><b>　　figure,</b></p><p>　　f=linspace(-fs/2,fs/2,nu

20、m); %頻率</p><p>　　plot(f,mag) %畫頻譜圖</p><p>　　title('Frequecy spectrum')</p><p>　　xlabel('Frequecy(Hz)')</p><p>　　ylabel('Magnitude(v/Hz)')<

21、/p><p>　　2.（1）利用MATLAB函數庫中函數設計butterworth高通濾波器，對采集的兩路信號分別進行高通濾波，抑制工頻噪聲（典型值50Hz）及其他常見低頻干擾噪聲。</p><p>　　MATLAB中的函數butter(N,Wn,’high’)實現(xiàn)butterworth高通濾波器，其中N表示濾波器的階數，Wn表示截止頻率。工頻噪聲為50Hz，由于巴特沃茲濾波器是非理想濾波器，

22、并且有其他頻率的低頻噪聲，所以截止頻率應高于50Hz，故這里截止頻率取100Hz。其中Wn應為0到1之間的數，Wn=截止頻率/采樣頻率的一半。通過filter這一維數字濾波器實現(xiàn)輸出。Y = filter(B,A,X) ，輸入X為濾波前序列，Y為濾波結果序列，B/A 提供濾波器的傳遞函數，B為分子， A為分母 整個濾波過程是通過下面差分方程實現(xiàn)的： a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... +

23、 b(nb+1)*x(n-nb) - a(2)*y(n-1) - ... - a(na+1)*y(n-na) [Y,Zf] = filter(B,A,X,Zi)，輸入X為濾波前序列，Y為濾波結果序列，B/A 提供濾波器系數，B為分子， A為分母。需要注意的是，這里的X和Y均為時域上的信號。用subplot使濾波前后的兩幅圖畫到同一張圖上。</p><p>　　濾波前后的頻譜圖如下：(上圖為原始信號，下圖為濾波

24、后的信號)</p><p>　　從頻譜圖對比，可以看出，信號的低頻分量已經被過濾掉了，濾波器設計成功。</p><p><b>　　具體代碼如下：</b></p><p>　　function Y=PreProcess(X) %X為原始時域信號</p><p>　　num=length(X);</p>

25、<p>　　[B,A]=butter(9,0.02,'high');%截止頻率100Hz</p><p><b>　　F=fft(X);</b></p><p>　　Y=filter(B,A,X); %其中X,Y均為時域信號</p><p><b>　　Z=fft(Y);</b></

26、p><p>　　Y1=fftshift(Z);</p><p>　　fs=1e4; %采樣頻率</p><p><b>　　figure,</b></p><p>　　f=linspace(-fs/2,fs/2,num);</p><p>　　subplot(2,1,1) %將濾波前后的信號畫到一個

27、圖上</p><p>　　plot(f,abs(fftshift(F))) %濾波前的頻譜圖</p><p>　　xlabel(‘Frequency(Hz)’)</p><p>　　subplot(2,1,2)</p><p>　　plot(f,abs(Y1)) %濾波后的頻譜圖</p><p>　　xlabel

28、(‘Frequency(Hz)’)</p><p>　　（2）將濾波后的兩組數據分別寫至data1p.txt和data2p.txt中。</p><p>　　MATLAB的文件操作跟C語言類似，用fopen打開一個文件，fprintf把數據寫入文件，fclose關閉文件。為了方便，這部分代碼寫在主程序下：（主程序的完整代碼見附錄）</p><p>　　%以下分別將濾波

29、后的兩組數據寫入data1p.txt 和data2p.txt</p><p>　　fid1=fopen('H:\學習\2012信號與系統(tǒng)-設計材料及要求\data1p.txt','w');</p><p>　　fprintf(fid1,'%12.8f\r\n',F1); %記事本中，只能用\r\n才能換行</p><p&

30、gt;　　fclose(fid1);</p><p>　　fid2=fopen('H:\學習\2012信號與系統(tǒng)-設計材料及要求\data2p.txt','w');</p><p>　　fprintf(fid2,'%12.8f\r\n',F2);</p><p>　　fclose(fid2);</p>&

31、lt;p><b>　　具體截圖如下：</b></p><p>　　用互相關方法檢測兩路信號的相對時間延遲，并以此計算漏點位置</p><p>　　互相關函數表示的是兩個時間序列之間和同一個時間序列在任意兩個不同時刻的取值之間的相關程度，即互相關函數是描述隨機信號x(t),y(t)在任意兩個不同時刻t1，t2的取值之間的相關程度。對于連續(xù)信號公式表示為</p

32、><p>　　R(τ)=（1/T）∫[f(t)g(t+τ)]dt，積分限為0至T。對于離散信號公式表示為</p><p>　　R(n)=（1/N）∑[x(m)y(m+n)]其中m從0到N-1變化。</p><p>　　MATLAB中，用互相關函數xcorr，畫出函數圖象。用max函數找出互相關函數的最大值所在的點。設原信號的長度為N，采樣頻率為fs，最大值點為Nmax，

33、則時間差D=(Nmax-N)/fs?？梢缘贸鲂孤c距離第一個傳感器的距離L1=（L-v*D）/2，到第二個泄漏點的距離為L2=(L+v*Dt)/2。</p><p>　　兩個信號y1和y2的互相關函數的圖象如下：</p><p><b>　　得</b></p><p>　　時間差 D=0.0121s</p><p>　　

34、泄漏點到第一個傳感器的距離 L1= 42.4375m</p><p>　　泄漏點到第二個傳感器的距離 L2= 57.5625m</p><p><b>　　具體代碼如下：</b></p><p>　　function Y=TDE(x1,x2)</p><p>　　Y=xcorr(x1,x2);</p>&l

35、t;p>　　Lx=length(x1);</p><p>　　Ly=length(Y);</p><p><b>　　fs=1e4;</b></p><p><b>　　figure</b></p><p>　　plot(Y), hold on</p><p>　　[Y

36、_max,N]=max(Y); %求出最大值及對應的橫坐標值</p><p>　　plot(N,Y_max,'ro') %標出最大值點</p><p><b>　　Dn=N-Lx</b></p><p>　　D=Dn/fs %延遲時間</p><p><b>　　L=100;<

37、/b></p><p><b>　　v=1250;</b></p><p>　　L1=(L-v*D)/2</p><p>　　L2=(L+v*D)/2</p><p><b>　　四、結論</b></p><p>　　利用MATLAB可以對信號進行時域和頻率域的分析，并

38、且畫出其時域圖象和頻譜圖。實際情況下找不到像理論課上所講的理想濾波器，因此，在實際情況下使用butterworth高通濾波器可以很好地過濾掉低頻噪聲，并且階數越高，越接近理想濾波器，濾波效果越好。需要注意的是filter函數中，輸入和輸出信號都是在時域上的，否則會出現(xiàn)問題。通過找出兩個信號的互相關函數的最大值點，可以求出兩個傳感器接收信號的時間差，從而求出泄漏點到兩個傳感器之間的距離。</p><p>　　通過這

39、次課程設計，我更加熟練地使用了MATLAB，對信號與系統(tǒng)的概念和物理意義有了更加深刻的理解，鍛煉了我分析問題和解決問題的能力。</p><p><b>　　五、參考文獻</b></p><p>　　[1].Alan V. Oppenheim, Alan S. Willsky, 劉樹棠譯.信號與系統(tǒng).西安交通大學出版社.2008.</p><p>

40、;　　[2].John R. Buck, Michael M. Daniel, Andrew C. Singer, 劉樹棠譯.信號與系統(tǒng)計算機練習——利用MATLAB.西安交通大學出版社.2006</p><p><b>　　六、附錄</b></p><p>　　由于MATLAB不能像C語言那樣，可以把多個函數寫到同一個程序里。故這里需要寫一個主程序。主程序代碼如下：

41、</p><p>　　%信號課程設計主程序</p><p><b>　　clear</b></p><p>　　y1=load('H:\學習\2012信號與系統(tǒng)-設計材料及要求\data1.txt');</p><p>　　y2=load('H:\學習\2012信號與系統(tǒng)-設計材料及要求\data

42、2.txt');</p><p>　　f=1e4; %采樣率</p><p>　　N=1/f; %采樣周期</p><p>　　FT1=SignalTF(y1,N);</p><p>　　F1=PreProcess(y1);</p><p>　　FT2=SignalTF(y2,N);</p>

43、<p>　　F2=PreProcess(y2);</p><p>　　%以下分別將濾波后的兩組數據寫入data1p.txt 和data2p.txt</p><p>　　fid1=fopen('H:\學習\2012信號與系統(tǒng)-設計材料及要求\data1p.txt','w');</p><p>　　fprintf(fid1,&#

44、39;%12.8f\r\n',F1); %記事本中，只能用\r\n才能換行</p><p>　　fclose(fid1);</p><p>　　fid2=fopen('H:\學習\2012信號與系統(tǒng)-設計材料及要求\data2p.txt','w');</p><p>　　fprintf(fid2,'%12.8f\r\