基于MATLAB的地震數字濾波器分析摘要對于地震而言，要降低財產損失和人員傷亡，最有效的方法是利用地震探測儀進行預警，地震探測儀的設計理念就是收集地震信號、轉換地震波、過濾地震波、分析地震波和發出預警。其中，過濾地震波是至關重要的一個環節，只要將地震波中的噪聲盡可能的過濾掉，留下的有效信號更準確，才更有利于后續的研究和分析。對于濾波器而言，對其性能和功能也提出了越來越高的要求，本課題的研究內容為基于MATLAB的地震數字濾波器。主要設計了Kalman和Wiener兩種濾波器，仿真模擬地震信號并加入一個15db和10db的白噪聲，比較濾波前后的頻率和幅值變化，通過計算信噪比得出Wiener濾波器比Kalman濾波器的降噪效果好。關鍵詞：地震信號;卡爾曼濾波器;維納濾波器;噪聲信號目錄TOC\o"1-2"\h\z\u引言 1第一章 緒論 11.1課題研究的內容及意義 11.2國內外研究的現狀 1第二章數字濾波器 22.1信號處理的發展和應用 22.2數字信號處理器 32.3數字濾波器 3第三章Wiener濾波器 43.1Wiener濾波器的簡介 43.2Wiener濾波器的應用 6第四章Kalman濾波器 74.1Kalman濾波器的基本理論 84.2Kalman濾波器的設計原理 10第五章MATLAB與濾波器的設計 115.1MATLAB的基本應用 115.2使用MATLAB語言進行程序設計 12參考文獻 18附錄A 19引言我國由于位于世界兩大地震帶之間，受三大板塊的時刻擠壓，所以是世界上地震災害最為嚴重的國家之一，其地震具有高頻率、高強度、淺震源等特點，這一列地震特點導致了地震災害的突發性、嚴重災害性[1]。但隨著電子工業的飛速發展，我國的地震探測器技術也在不斷增強，但對于地震探測儀而言準確的地震波信號是至關重要的一部分。本次實驗主要針對地震信號的數字濾波器進行研究，使用MATLAB仿真軟件進行仿真，使用卡爾曼和維納兩種濾波方法設計濾波器，輸入模擬地震信號并加入不同頻率的白噪聲，對比濾波前后幅值和頻率波形的變化，計算信噪比，最后通過比較信噪比大小得出維納濾波器的效果更好。緒論1.1課題研究的內容及意義隨著技術的不斷發展，數字信號處理在許多領域幾乎已經取代模擬信號，人們對它轉換的精度的要求也越來越高。電子信息類的工業發展非常迅速，對于濾波器而言，它的功能和性能的要求也越來越多。雖然在以前傳統的地震勘探儀器中模擬濾波器使用非常廣泛，但是他也存在很多問題，一方面在使用過程中會出來偏移現象，另一方面還會存在器件噪聲等等，這些缺點在一定程度上影響了對地震勘探的工作。但是現在數字信號處理技術正在飛速發展，數字濾波器以它穩定的性能、簡單的操作、傳輸方便等特點淘汰了傳統的模擬濾波器。1.2國內外研究的現狀隨著國內經濟的不斷發展，祖國在建設方面的發展速度也越來越快。但地震對人類的威脅還是不可忽略的。為了將地震引起的傷亡和損失降到最低，地震勘探儀器起著不可忽略的作用。隨著科技的進步，地震探測儀也經過了長達二十年的發展，已經實現了數字勘探設備。對于地震勘探儀而言數據采集的數字部分是不可或缺的一個模塊，只有將地震信號更好的檢測并進行濾波，得到去除掉干擾波的有效波，才能使得地震勘探更加精確。地震信號是一種振動信號，他有兩個特點，一個是它的超過120db的動態范圍，另一個是高分辨率。地震勘探主要分為三個階段：數據收集、數據處理和地震材料分析[2]。我們可以通過處理地震數據來獲得與地震有關的信息。如何去除信號中的干擾噪聲，得到有用的地震波形數據是信號處理過程中的重點、難點[3]。目前，對于地震的勘探而言，設計勘探儀器最主流的方法是利用地震法，地震法的原理是利用檢波器來采集到所需要檢測的地震波，地震波并不會產生后直接被收集起來，而是在地層中首先要經過傳播，其次是折射，最后是反射，然后再返回到地表面，之后再對這些數據進行研究和分析。通過以上介紹可知，當前國際主流的地震信號濾波器方案已經很成熟，性能優秀，廣泛應用于萬道量級及以下的勘探儀器中，能夠滿足采集需求[4]。當前地震勘探儀器由負責數據采集的采集站、負責電源供給和數據轉發的電源站、負責數據匯總的交叉站和負責記錄數據和控制站體的主機構成[5]。對于采集站而言，數字濾波器對于它的設計具有巨大的作用。第二章數字濾波器2.1信號處理的發展和應用在大自然中，信號無處不在，信號的種類有很多，可以分為自然產生的和有意制造出來的，也可以分為有用的和沒用的等，無論是哪一種信號，人們通過數字信號處理都可以提取出有用的信號。隨著電子工業和信息工業的不斷發展，數字信號處理技術已經成為了一個極其重要的問題和技術。數字信號對于模擬信號而言具有以下幾點優點。靈活性大。在進行程序的編寫和修改時，數字濾波器只需要改變部分程序、代碼或者數字就可以實現用戶的需求。精度高。數字模電路的元件精度要比模擬電路的高很多。可靠性高。數字系統有且只有“0”、“1”兩個電平信號，不易被外界影響而改變。采用大規模集成電路，所以其可靠性也相對較高。易于大規模集成。數字濾波電路的工作主要在截至或飽和狀態。并行處理等優點。所以在越來越多的領域中，數字信號發展的十分迅速，并且已經慢慢取代了模擬信號。2.2數字信號處理器數字信號處理就是對信號進行檢測和濾波，其算法不僅需要利用數字信號處理器，還需要使用不止一個專用集成電路進行處理。數字信號處理器如同數字信號一般，也具備靈活性高、精確性好等一系列的優點。其在信號處理的過程中一般需要進行三個步驟：信號的過濾與檢測、參數的提取和估計、頻譜分析等，通過這三個步驟才能使信號變得更加有利于人們的識別。2.3數字濾波器在數字信號處理中，濾波器占有極其重要的地位，數字濾波是語音和圖像處理、模式識別、譜分析等應用中的一個基本的處理算法[6]。在各個方向的研究中，人們不單單看重信號處理的快速性，同樣對實時性也有很大的要求。在對信號進行處理時，無論是過濾、檢測還是預測的過程中，都離不開濾波器的使用，數字濾波器的使用意義無非就是一種離散的時間系統可以把有限精度算法實現。數字濾波器是采用有限精度算法實現離散時間系統，通常應用于修改或者改變頻域中信號的屬性[7]。數字濾波器具有以下的優點：穩定性高。精度高。設計靈活。實現方便。其中數字濾波器具有穩定性高、精度高、設計靈活、實現方便等許多突出的優點，避免了模擬濾波器所無法克服的電壓漂移、溫度漂移和噪聲等問題，因而隨著數字技術的發展，用數字技術實現濾波器的功能越來越受到人們的注意和廣泛的應用[8]。應用最多的數字濾波器是線性時不變（LinearTime-Invariant，LTI）濾波器，LTI濾波器通常分為有限脈沖響應（FiniteImpulseResponse，FIR）濾波器和無限脈沖響應（InfiniteImpulseResponse，IIR）濾波器[9]。隨機信號和隨機過程的存在是非常普遍的。首先，每一個確定性信號在經過采集和測量的過程中總會不可避免地引入一些隨機信號，這種隨機信號可以被看作是一種隨機誤差，但是如果不存在隨機誤差，那么也就無法將確定性隨機化了。其次，每一種信號都無法避免隨機信號的干擾，所以給這種干擾賦予一個非常形象的稱呼——噪聲。但是對于干擾和噪聲我們又有兩種不同的概念，我們可以將干擾看作是一種確定信號也可以看作是噪聲，也就是隨機信號，它的組成可以分為兩個部分，一部分是白噪聲一部分是確定性信號。那么區分兩種信號通常是判斷是否為目標信號，如果是目標信號則為噪聲信號，如果不是目標信號，那么就是干擾信號。因此對于濾波器而言，一般是由確定信號伴隨這干擾信號同時輸入，并且輸出時應該盡可能的去除掉干擾信號。維納濾波和卡爾曼濾波就是解決這類問題的兩種方法[10]。噪聲可以按照功率譜密度進行劃分可以分為以下兩種。白噪聲（WhiteNoise），他也被稱作隨機信號。色噪聲（ColorNoise）。每一種隨機信號都可以看作是由白噪聲和確定性信號隨機組合的。最優濾波主要解決系統的狀態或信號的最優估計問題,即由被噪聲污染的觀測信號求在某種性能指標和某種意義下狀態或信號的最優估值器，也叫最優濾波器[11]。最優的定義是相較而言的，一般情況下對于濾波器選用線性濾波器，因為輸入的信號通常是線性函數，對于濾波的定義無非是過濾波形，因為在測量信號時會有一些噪聲的存在，噪聲會影響對于波形的判斷，除了噪聲儀器的使用也會導致出現一些誤差影響信號的測量結果。通常解決濾波的三種方法分別是：Wiener濾波方法。使用維納濾波的基本方法可以分為兩種一種是利用譜分解的經典維納濾波，另一種是利用Diophantine方程的現代維納濾波。Kalman濾波方法。基本方法是利用Riccati方程。現代時間序列分析方法應用現代時間序列分析方法，基于ARMA新息模型和白噪聲估值器[12]。本次實驗主要使用了前兩種濾波方法分別對模擬地震波進行濾波，并選取最適合地震波的濾波方法。第三章Wiener濾波器3.1Wiener濾波器的簡介維納濾波方法是由N.Wiener在20世紀40年代初由于研究火炮控制系統的需要提出來的[13]。它使用的是平穩隨機過程譜分析，而信號的平穩隨機又成為了它的局限，所以不得不存儲所有的歷史數據，對于非遞推濾波器而言，他不僅存儲量非常的巨大，它的計算量也是不容小覷的，所以以上種種缺點導致了他存在不方便使用的問題。但是就算經典維納濾波由上述種種缺點和問題，但它對于濾波器的設計依舊是非常重要的工具和方法。上世紀八十年代后，又出現了現代維納濾波方法，他對于經典維納濾波而言是一種新的突破，它不再局限于平穩隨機的信號，而是可以處理非平穩非隨機，不僅如此，它還可以處理多維信號，它突破了經典維納濾波方法的局限。維納濾波的重要特性有：正交性原理。誤差信號統計表征為白噪聲[14]。假設，存在一個線性系統，h(n)代表單位脈沖響應，x(n)為輸入的一個隨機信號，它包含噪聲w(n)和有用信號s(n)，則：X(n)=s(n)+w(n)(3-1)則輸出為：yn=xn為了使輸出信號最大程度接近有用信號，稱輸出信號為有用信號的估計值s(yn=syx(n)=s(n)=w(n)h(n)圖3-1維納濾波器的系統框圖濾波的定義就是用當前的和過去的觀測值來估計當前的信號：yn=s(預測，顧名思義就是用過去的觀測值來估計當前的或將來的信號：yn=s(平滑或者內插，簡單點說就是用之前的觀測值來推測過去的信號：yn=sn+上圖中的估計信號和想要得到的有用信號是不可能完全相同的，在這里用e(n)來表示真值和估計值之間的誤差：en=sn雖然e(n)是隨機變量，維納濾波和卡爾曼濾波的誤差準則就是最小均方誤差準則：E[e2(n3.2Wiener濾波器的應用設計維納濾波器需要做幾個方面的準備，首先需要估計統計特性，然后再了解清楚觀測信號和估計信號直接的關系，以及用到的一些函數，最后才能設計出維納濾波器，這種設計方法設計出來的濾波器又叫做——后驗維納濾波器。維納濾波器應用領域十分廣泛，例如在生物醫學信號處理中非常常見的應用就是有關提取誘發腦電信號的相關研究。研究者在提高信噪比時采用了維納濾波的濾波方法，前后有多名科學家對濾波方法進行了研究和改進。對于頻域應用后驗維納濾波而言，他最重要的地方就是需要使用每一次觀測信號分解出來的信號和噪聲的譜估計，最后可以通過使用設計的濾波器提高信噪比。在白噪聲情況下最有濾波算法本質上是一個利用最小二乘法進行最小方差擬合的過程[15]。考慮真實世界里的儀器，對真實信號的響應(response)不是完美的delta函數,而是由真實信號u(t)與某個已知的響應函數r(t)卷積所得的“被彌散化”(smeared)的信號s(t)：st=?∞∞寫成傅立葉變換形式為：S(f)=R(f)U(f)(3-9)另外,儀器所測得的信號c(t)還含有噪音成分n(t)，因此c(t)=s(t)+n(t)。c(t)也稱為損壞(corrupted)的信號。在沒有n(t)時,我們可以使用去卷積算法,利用傅立葉變換，能很容易地計算得到真實信號：U(f)=S(f)/R(f)(3-10)而此時需要考慮有n(t)時，如何得到u(t)。此時，設存在一個最優濾波器?(t)orΦ(fUf=C(f)Φ(f)R為了盡量接近真實信號，則使差的模平方最小即可：st=?∞∞將（1）式代入（2）式，得到：?∞∞|Sf對上式，考慮到模的平方，有個小技巧：因為信號與噪音無關(uncorrelated)，所以含有交叉項S(f)*N(f)的積分為0.將S(f)和N(f)分離開，展開模平方，得到：?∞∞1Rf我們要找到一個對測量值有意義的實函數Φ(f),使上式的積分核最小化。因此，將其他量視為常數，對Φ(f)求導并令結果=0,同時考慮到Φf(1?Φf)2'=2可以求得：Φf=|S(f此即最優濾波器，或者Wiener濾波器的公式。這個公式里不含U和R，只需要估計S^2和N^2。當噪音可以忽略時Φf→1,當噪音占主導時，第四章Kalman濾波器維納濾波理論是四十年代美國科學家Wiener等人發現的，它又有另外一個名稱是最佳線性濾波理論。一般來說，維納濾波在使用過程中存在的最大的難題就是需要使用之前所有的數據，但是這個要求應用在需要隨時處理的信號而言是很困難的。為了解決實時性的問題，卡爾曼在六十年代的時候將狀態空間模型反映到濾波理論中，并設計出一套遞推估計算法，也就是我們所說的卡爾曼濾波理論。基本思想是：采用信號和噪聲的狀態空間模型，利用前一時刻的估計值和現時刻的觀測值來更新對狀態變量的估計，求出現時刻的估計值[16]。它的出現不僅彌補了維納濾波的缺點，而且應用范圍也更加廣泛了，可以應用在實施運算和計算機運算中。隨著電子工業的迅速發展，經典維納濾波方法已經越來越無法滿足實際應用的需求。對于濾波算法的要求最好是遞推的，這樣方便進行實時計算，然后可以處理多變量非平穩隨過程或時變系統濾波問題。卡爾曼濾波是一種時域濾波方法，它解決最優濾波問題的理論基礎是狀態空間模型和射影，然后再通過計算或求解Riccati方程。卡爾曼濾波的優點是濾波算法可遞推，這樣方便在計算機上實現，它應用在多個場合中，例如時變系統、非平穩信號和多維信號濾波問題等，現代維納濾波最大程度解決了經典維納濾波的局限性和其他缺點，但它并不是完美的，它也需要提前精確了解系統模型和噪聲方差。4.1Kalman濾波器的基本理論卡爾曼（Kalman）濾波，是一種非常常見的濾波算法，相對于Wiener濾波而言它具有很多的優點。通過不斷的預測和實測來修正自己的估計值，最后達到一個理想的平穩狀態[17]。狀態空間方法是卡爾曼濾波方法的基本特征和關鍵技術，狀態空間模型是卡爾曼濾波器的設計基礎，而狀態空間模型的基本概念是系統的狀態變量，系統狀態變量是描述系統特征的n維列向量，它的取值空間為n維歐氏空間。例如，假設有一個做直線變速運動的點，它的狀態通常使用位置、速度和加速度三個要素來表示。如果狀態與加速度無關，則可以省略加速度這一要素。狀態變量對比信號而言，它的概念更廣泛。根據不同的情況有不同的定義，比如說可以把信號看成狀態或者狀態的分量。狀態空間模型一方面是狀態空間方法基本模型，另一方面是出發點。狀態空間模型又包括兩種方程其一是狀態方程，其二是觀測方程。狀態方程又稱作狀態模型，他通常用來描述狀態變化的規律；而觀測方程又稱作觀測模型，它用來定義對狀態進行線性觀測的方程，需要注意的是觀測方程一般情況是含有隨機觀測噪聲。對于線性觀測的定義是被觀測信號可以是狀態本身，也可以是狀態的一種或多種分量，也可以是狀態分量的線性組合。但觀測方程也可以非線性的。一般噪聲都是在全頻段都有分布的，這是與有用信號處于同一頻帶的噪聲就無法通過FFT算法過濾了[18]。在卡爾曼濾波算法中,首先對某時刻系統的擾動和觀測誤差(即噪聲)展開統計,然后在線性狀態空間表示的基礎上,通過處理含噪的觀測信號求得誤差最小時真實信號的估計值,從而排除數據中噪聲和干擾信息的影響[19]。卡爾曼濾波適用于所有含有不確定因素的動態系統，可以通過使用一些數學建模對系統接下里的動作作出大致的預測，雖然很難保證系統不會受到外界的干擾，但是卡爾曼濾波在絕大多數情況下都可以提高系統預估的精確度，然后就可以清晰的知道系統狀態的轉移情況。并且卡爾曼濾波還可以有效的利用多個粗糙數據之間的關系。卡爾曼濾波尤其適用于動態系統，因為它存在以下兩個優點：它對于內存的要求很低，因為它僅僅需要保留系統上一個狀態的數據，而不需要保留之前很久一段的數據。它的運算速度很快，非常適合解決實時問題和應用于嵌入式系統中。卡爾曼濾波其實是這樣一個線性高斯系統應用貝葉斯濾波方程而得到的顯式解。hiddenxR+1xRxR-1hiddenxR+1xRxR-1YR+1YRYR-1observedYR+1YRYR-1observed圖4-1卡爾曼顯示解圖形卡爾曼濾波在進行濾波時有兩個步驟：第一步是預測，系統根據前一個狀態的結果對下一個狀態進行判斷和估計，第二步是更新，根據目前的結果實時更新系統的狀態。對于Kalman而言有五個公式是非常常用的。1、狀態先驗估計xk|k?1這個公式就是我們研究的變量在這一時刻與上一時刻的具體關系。從這個公式也可以看出卡爾曼濾波只能應用于線性系統中的。2、方差先驗估計Pk|k這個公式是用來計算估計值與實際值的方差大小，通常可以利用這個公式評估當前估計的誤差有沒有問題。3、增益矩陣Kk=Pk卡爾曼濾波通過估計當前觀測變量的變化趨勢來估計某一時刻變量的值，這個增益矩陣就是指規定的這個變化量的大小。4、狀態估計校正xk|這個式子是用來求估計的值。5、方差估計矯正Pk|k=[這個公式是用來修正估計值和實際值的方差。卡爾曼濾波解決的是如何從多個不確定數據中提取相對精確的數據。實驗前提是這些數據必須滿足高斯分布，而理論前提是第一個高斯斑乘以第二個高斯斑從而得到第三個高斯斑，而這最后一個高斯斑就是提取到相對精確的一個數據范圍。卡爾曼濾波器（KalmanFilter,KF）類算法是由美國學者Kalman在20世紀60年代初提出的，此類算法是一種最小方差意義上的最優估計方法[20]。與此同時，卡爾曼濾波器也是一種有關最優化回歸數據的處理算法。在很多問題的處理方面，它可以算得上是最優的并且效率也是最高的，毫不夸張的講它是非常有用的。在過去三十多年的時間段里，它的適用領域不僅僅在日常生活中（人臉識別），還存在于導航領域和控制領域，甚至在軍事領域（雷達）都被廣泛應用。目前卡爾曼濾波算法在各種復雜情況和環境下的目標跟蹤問題上任然是研究和應用熱點[21]。4.2Kalman濾波器的設計原理卡爾曼濾波器使用狀態觀測器。狀態觀測器是針對可觀系統，根據輸出y和輸入u對系統內部狀態u進行觀測的結構單元。其構圖如下：x=Ax+Buxy=Cx+_+_x=Axxyyxxkk圖4-2狀態觀測器結構圖原系統方程為：

x=Ax+Buy=Cx(4-7)采用狀態觀測器的觀測系統方程為：x=Ax+y=Cx為保證觀測器的limt→∞x?卡爾曼濾波器就是一種狀態觀測器，只不過它是隨即系統的狀態觀測器，其結構框圖如下：KalmanfilterKalmanfilteryvyyeyvyyeuplantuplantSensornoiseProcessnoiseySensornoiseProcessnoisey圖4-3卡爾曼濾波器原理的結構框圖在輸入u和動態系統plant中間會引入過程噪聲w，而在輸出y和實際測量yv中間會引入傳感器噪聲v，而卡爾曼濾波器則是根據u，yv求得測量的最優估計ye。第五章MATLAB與濾波器的設計5.1MATLAB的基本應用5.1.1MATLAB的簡介MATLAB是一個高效率的數值計算可視化軟件，它的全程叫做“矩形實驗室”(MatrixLaboratary)，是美國Mathworks公司開發的一種大型數字計算的軟件，他不僅應用范圍非常廣泛，包括工業、電子、醫療、建筑等領域，例如：信號和圖像的處理、檢測和分析等，他的優點也十分顯著，包括開放性優良、靈活性高、使用簡單等。MATLAB在世界范圍內也成為了一種非常流行的軟件，被越來越多的工作人員、老師和使用。MATLAB語言一方面可以被看成是解釋性語言，因為它擁有十分強大的語言編程和交互式計算機的環境，操作者不僅可以利用軟件發出一些命令使程序進行操作或數字運算，而且也可以通過使用不同的語言對程序進行編寫、運行和調試。另一方面MATLAB語言是一種開放性的語言，用戶可以通過實時編寫和修改程序以此達到自己的目的。5.1.2MATLAB的特點1、MATLAB在科學計算方面具有完善的工具，所有的變量都是矩陣對象，因為采用矩陣運算而不是循環運算，所以速度快。2、最接近通用語言的科學計算語言。3、支持各種語言的擴展，例如python、c語言等等。4、語法簡單，便于操作。5、并行實現較容易。5.2使用MATLAB語言進行程序設計使用模擬地震波，在模擬地震波的基礎上加入不同分貝的白噪聲，再使用兩種濾波器對其濾波，分別對比濾波后和原模擬地震波的時間-幅值和頻率-幅值波形的變化，計算出信噪比并進行比較分析，選擇效果最好的濾波器。（程序見附錄A）Matlab的程序主要分為幅值和頻率兩個部分，在編寫程序時分為以下幾個步驟：讀取地震波。使用dire函數指定路徑。使用filename函數將變量導入工作區，以便后續的使用。使用stream隨機流將模擬地震波的數據與時間一一對應畫出時間-幅值函數和頻率-幅值函數。使用data做出矩陣。對函數做傅里葉變換（fft），生成時間-幅值和頻率-幅值圖。加入隨機白噪聲。隨機設置白噪聲的SNR，單位為db。使用randn(size(data1))函數。利用mean函數求噪聲平均值，再用原始噪聲減去其平均值，得到噪聲信號。計算出信號的的強度后計算出噪聲的方差。利用噪聲方差的平方根/信號的標準差乘以信號生成噪聲信號。利用awgn函數，y=awgn(x,snr,'measured')和y=awgn(x,snr)是相同的，在添加噪聲之前測量了x的能量。給正弦波信號加入高斯白噪聲生成觀測信號。定義誤差信號為加入的白噪聲（也就相當于觀測信號減去地震模擬信號。生成圖像。使用卡爾曼濾波算法進行濾波。測量更新預測更新預測誤差協方差矩陣測量更新預測更新預測誤差協方差矩陣預測誤差協方差矩陣預測誤差協方差矩陣狀態量估計狀態量估計預測狀態量更新K增益矩陣預測狀態量更新K增益矩陣圖5-1卡爾曼濾波算法實現流程圖設置初始化變量后，計算并賦值轉移和測量矩陣、系統噪聲矩陣、測量噪聲