1、分類號 密級本 科 畢 業 設 計題 目經驗模態分解在地球物理資料中的應用 英文題目 Empirical Mode Decomposition and its Application in Geophysics 二OO七年六月本科畢業設計論文任務書學生姓名班級專業電子信息工程導師姓名職稱單位畢業設計論文題目經驗模態分解在地球物理資料中的應用畢業設計論文主要內容和要求：主要內容：1、研究經驗模態分解法原理；2、利用MATLAB語言編程實現經驗模態分解法；3、對理論模型和實際數據進行試驗，驗證經驗模態分解法。要求：經驗模態分解法是一種新的分析非線性非平穩信號處理方法。利用經驗模態分解法處理地球物理

2、非線性非平穩信號，別離提取資料中的特征信息，進而可對其分析解釋。畢業設計論文主要參考資料：1、Hassan H. H., Empirical Mode Decomposition (EMD) of potential field data, SEG 72th Annual International Meeting Expanded Abstracts, 2005: 704-706.2、曾華霖，重力場與重力勘探，地質出版社，20053、管志寧，地磁場與磁力勘探，地質出版社，2005畢業設計論文應完成的主要工作：1、研究經驗模態分解法原理；2、用MATLAB語言編寫經驗模態分解法程序；3、理論模

3、型和實際地球物理資料的經驗模態分解法試驗；4、分析總結數據試驗結果。畢業設計論文進度安排：序號畢業設計論文各階段內容時間安排備注1課題調研、查詢資料、外文資料翻譯、文獻綜述，準備開題報告3月5日3月30日2經驗模態分解法方法研究和程序編寫實現，準備中期檢查4月1日5月7日3理論模型和實際數據試驗，改良完善方法和程序5月8日5月20日4撰寫論文，準備辯論5月21日6月3日課題信息：課題性質： 設計 論文 課題來源： 教學 科研 生產 其它發出任務書日期： 2007年3月10日 指導教師簽名： 年 月 日教研室意見：教研室主任簽名： 年 月 日學生簽名： 摘 要經驗模態分解(EMD)是由Huang

4、等人提出的一種新的分析非線性、非平穩信號的方法。本文研究經驗模態分解原理及其在地球物理資料中的應用。首先研究經驗模態分解的根本原理和算法，對地球物理資料地震資料，重磁資料進行EMD分解試驗分析，然后研究基于EMD的Hilbert變換原理及其在提取地震屬性信息中的應用，對實際地震時間剖面和時間切片進行EMD時頻分析試驗。本文的方法研究和數據試驗分析說明：經EMD分解變換得到的IMF序列是直接從原始時序數據中別離出來的,事先無需確定分解階次,能更好反映原始數據固有的物理特性,每階IMF序列都代表了某種特定意義的頻帶信息；EMD分解獲得的IMF序列具有穩態性,對IMF進行Hilbert變換，就可以得

5、到單個固有模態函數的瞬時振幅、瞬時相位和瞬時頻率，這些信息可以清楚的顯示信號的時頻特征；EMD分析方法用于分解地球物理資料和作時頻分析是有效的。關鍵詞：經驗模態分解；地球物理；Hilbert變換；固有模態函數；時頻分析ABSTRACTEmpirical Mode Decomposition(EMD), which was developed by huang, is a new method to analyse nonlinear and nonstationary signals. In this paper, we study the theory of EMD and its appl

6、ications in handling geophysical data. Firstly, we introduce the theory and the Methodology about EMD ,then we will use this method to analyse the geophysical information, including the gravity anomaly data and seisms data. Based on the EMD, we will study the theory of the Hilbert transform, and the

7、n use it to obtain the images,from which we can deal with the seisms slice by time- frequency analysis in order to distill the seisms information. The studying of EMD and the data testing in this paper indicate: intrinsic mode functions(IMF) is comes from the original signal by the EMD, in this cour

8、se, we need not fix on the Decomposition number and would not influenced by some mens factors. Every intrinsic mode function stand for some given information and can reflect the intrinsic physical Information very well. Meantime, take the IMF for Hilbert transform, then we can get the IMF frequency

9、and the amplitude, which can show the signals characteristics very clearly. So it is very useful to use the EMD to study the geophysical data.Key words: Empirical Mode Decomposition; geophysics; Hilbert transform; intrinsic mode function; time- frequency analysis目 錄 TOC o 1-3 h z u HYPERLINK l _Toc1

10、69234403 第一章 緒 論 PAGEREF _Toc169234403 h 1 HYPERLINK l _Toc169234404 第二章 經驗模態分解(EMD) PAGEREF _Toc169234404 h 2 HYPERLINK l _Toc169234405 第一節 EMD根本原理 PAGEREF _Toc169234405 h 2 HYPERLINK l _Toc169234406 第二節 EMD算法流程 PAGEREF _Toc169234406 h 3 HYPERLINK l _Toc169234407 第三節 EMD例如 PAGEREF _Toc169234407 h 5

11、 HYPERLINK l _Toc169234408 第三章 EMD分解地球物理信號 PAGEREF _Toc169234408 h 8 HYPERLINK l _Toc169234409 第一節 EMD分解位場資料 PAGEREF _Toc169234409 h 8 HYPERLINK l _Toc169234410 一、EMD分解一維重力異常數據 PAGEREF _Toc169234410 h 8 HYPERLINK l _Toc169234411 二、EMD分解二維重力異常數據 PAGEREF _Toc169234411 h 9 HYPERLINK l _Toc169234412 第二節

12、 EMD分解地震資料 PAGEREF _Toc169234412 h 13 HYPERLINK l _Toc169234413 一、EMD分解一維地震數據 PAGEREF _Toc169234413 h 13 HYPERLINK l _Toc169234414 二、EMD分解地震時間剖面 PAGEREF _Toc169234414 h 14 HYPERLINK l _Toc169234415 三、EMD分解地震時間切片 PAGEREF _Toc169234415 h 17 HYPERLINK l _Toc169234416 第三章 基于EMD的Hilbert變換與地震信號屬性信息 PAGERE

13、F _Toc169234416 h 19 HYPERLINK l _Toc169234417 第一節 基于EMD的Hilbert變換原理 PAGEREF _Toc169234417 h 19 HYPERLINK l _Toc169234418 第二節 提取地震信號屬性信息 PAGEREF _Toc169234418 h 20 HYPERLINK l _Toc169234419 一、提取地震時間剖面屬性信息 PAGEREF _Toc169234419 h 20 HYPERLINK l _Toc169234420 二、提取地震時間切片屬性信息 PAGEREF _Toc169234420 h 22

14、HYPERLINK l _Toc169234421 結 論 PAGEREF _Toc169234421 h 23 HYPERLINK l _Toc169234422 致 謝 PAGEREF _Toc169234422 h 24 HYPERLINK l _Toc169234423 參考文獻 PAGEREF _Toc169234423 h 25 HYPERLINK l _Toc169234424 附錄：經驗模態分解MATLAB程序 PAGEREF _Toc169234424 h 26第一章 緒 論目前，信號處理1方法主要有傅立葉變換、小波變換等，這些方法在處理線性、平穩的信號時，具有簡單、效率高等

15、優點，然而，在實際的地球物理數據測量中，許多信號是非平穩的或非線性的，如果在處理這些信號時，仍假設數據是平穩或線性的，會使我們得到錯誤的分析結果，為了正確有效地分解非平穩信號，本文研究一種新的方法，即經驗模態分解(EMD) 2，該方法是由Huang等人提出的一種新的分析非線性、非平穩信號的方法。經驗模態分解EMD方法可將信號分解為假設干個固有模函數(Instrinsic Mode Function，IMF)，這些IMF是完備的并且正交的。由于EMD方法是依據數據本身的時域信息進行的時域分解，得到的IMF的個數通常是有限的和平穩的，基于這些IMF分量進行的Hilbert變換，得到IMF的瞬時振幅

16、，瞬時頻率，以及瞬時相位，從這些頻譜圖上我們可以得到真實清楚的物理信息，因此，其Hilbert譜也能夠準確反映出信號能量、頻率在空間或時間尺度上的分布，在非線性和非平穩過程的分析中具有很高的應用價值。由于經驗模態分解法主要用于分析非線性非穩定性數據，所以它在各種資料處理中應用非常廣泛。在過去的幾年,EMD方法已經被廣泛用于非線性的海洋波動數據分析3、地震信號和結構分析4、位場數據分析5、機械故障診斷6、醫學信號7、橋梁和建筑物狀況監測的分析、太陽輻射的變化分析等領域的研究。也可以將經驗模態分解法擴展至二維，應用于圖象處理8，將圖像分解為一系列不同尺度的信息，為深入處理圖象提供足夠的不同尺度的信

17、息，在實際應用中取得了很好的效果。本文研究經驗模態分解原理及其在地球物理資料中的應用。本文首先研究經驗模態分解的根本原理和算法，對地球物理資料地震資料，重磁資料進行EMD分解試驗分析，然后研究基于EMD的Hilbert變換原理及其在提取地震屬性信息中的應用，對實際地震時間剖面和時間切片進行試驗分析。第二章 經驗模態分解(EMD)第一節 EMD根本原理經驗模態分解法利用信號的局部特征時間尺度，從原信號中提取出假設干個固有模態函數(IMF)和一個剩余量，分解出的各個IMF分量突出了數據的局部特征，剩余分量表達了信號中的緩慢變化量，對它們進行分析，可以更準確有效地把握原數據的特征信息。用EMD來分解

18、一個非平穩非線性的信號是建立在一個假設根底之上的:任意信號都是由有限個固有模態函數(IMF)組成的。而固有模態函數具有兩個特征:(1)整個信號中的過零點和極值點的個數相等或相差1;(2)信號關于時間軸局部對稱，即信號任意一點由局部極大值確定的包絡和由極小值確定的包絡的均值為零。EMD方法將信號分解成假設干個IMF的過程如下：先找出原始信號(t)的局部最大值或最小值，用三次樣條函數進行插值，分別形成上包絡和下包絡。然后再求上、下包絡的平均值和差值，即：上下包絡的平均值公式如下： 1-1原函數(t)和的差值公式如下 (1-2)如果中極值點的數目和跨零點的數目相等或至多只差一個,并且各個瞬時平均值都

19、等于零,那它就是固有模函數，設IMF分量為，那么有 (1-3)否那么,把當作原序列,重復以上步驟,直至滿足固有模函數的定義,求出為止。第一個IMF分量代表原始數據中最高頻率的分量。將原始數據序列(t)減去第一個IMF分量可以得到一個去掉高頻分量的差值數據序列。對再進行上述平穩化處理過程可以得到第二個IMF分量,如此重復下去直到最后一個差值序列不可分解為止，此時代表原始數據的趨勢或均值：整個EMD過程用公式表示如下： 1-4 。 1-5原始數據序列可由這些IMF分量和一個剩余量疊加重構，公式如下： 1-6由于每一個IMF分量都是代表一組特征尺度的數據序列，因此這個平穩化的處理過程實際上是將原始數

20、據序列分解為各種不同特征波動的疊加。每一個IMF分量既可以是線性的也可以是非線性的。這種算法其實可以看作為一個篩選的過程，首先將信號中的頻率最高的成分篩選出來，在從原信號中將該成分去除，再從新的信號中選出頻率最高的成分，依次類推，直到信號不可分解為止。這種過程可以看作為一系列的濾波器族.從算法的以上描述可知，信號可以通過有限的篩選步驟完成分解。因為每次篩選時，新信號的最大、最小值的個數都在減少。為了減少提取IMF的篩選步驟，定義了參數: 1-7當小于某一常數時停止篩選，一般SD的值在至之間。另外在篩選過程中，由于該算法采用的是三次樣條插值，所以當信號的極大值或極小值的個數小于2時，停止篩選。這

21、些有限的IMF經過希爾波特變換后產生了以時間為參數的具有實際意義的瞬時頻率。這樣可以在時域和頻域同時對信號進行有效的分析，第二節 EMD算法流程根據對EMD過程的論述，EMD的詳細的算法如下：1.初始化：，i=1;2.得到第i個IMF:(a)初始化: =1,(b)找出的局部極值點，(c)對的極大和極小值點分別進行插值，形成上下包絡線,(d)計算上下包絡線的平均值(e)(f)假設0.3，那么，否那么j=j+1,轉到b3.的極值點數不少于2，那么i=i+1,轉到2，否那么分解結束，是剩余分量。EMD算法的流程圖如下：YESYESNO開始C=Xh=C,SD=1找出h的局部極值點極值點數2嗎？對局部極

22、大值極小值分別進行三次樣條插值，得到上下包絡線后求的包絡平均值mPrevh=h;h=h-m;根據公式計算SDSD0.3?保存當前的h為固有模函數H的極值點數2嗎？結束C=C-hNONOYES第三節 EMD例如本節以一單道地震數據為例，用于說明EMD方法處理數據的過程，以及EMD方法在處理非平穩非線性信號時所具有的特點。圖1為某非平穩非線性性一維地震數據圖。圖1 單道地震數據圖2顯示了圖1中信號的上包絡以及下包絡，m(t)為上下包絡的平均值。圖2 原信號的上下包絡以及包絡均值從原信號中分解IMF1的過程如下：圖3 從原函數中提取IMF1按以上方法，信號經EMD處理的結果如下：圖4 地震波信號的E

23、MD處理過程從圖4中，我們可以看出信號經EMD程序的平穩處理過程的特點：1頻率最高的成分最先被分解出來，IMF8缺乏一個周期，為剩余分量。2在振幅上，各個固有模函數也表達了明顯的變化：一般情況下，隨著信號EMD分解過程的進行，后面的固有模函數的幅值要逐漸變小。在圖四中，從IMF2到IMF8清楚的反映了這一變化。3圖4中是全部的固有模函數的疊加，通過與原函數相比，兩者根本上完全一樣。這也說明，運用EMD方法處理信號的整個過程不會造成信號的失真。4圖中整個EMD處理過程說明：我們可以通過疊加不同的固有模函數，對不同的固有模函數在原信號中所起的作用進行分析，也可以對某一固有模函數進行分析。比方在圖4

24、中，IMF5到IMF8，這幾個成分幅值小，頻率低，在原信號中并不占主導作用，這些成分也許是真實的某個過程中的現象，或是原數據中的噪聲，通過對其他固有模函數進行加權，就能濾除這些噪聲。需要說明的是，EMD方法是依據原函數本身的時域特征對原信號進行分解的，因此EMD方法具有有效性，真實性等特點，這也說明了EMD是一種非常好的非線性非穩定性的信號分析方法。第三章 EMD分解地球物理信號第一節 EMD分解位場資料在重磁勘探中，位場數據的有效性要受到數據獲得過程中所受噪聲干擾程度的限制，因為數據是非線性非穩定性的，所以我們采用經驗摸態分解EMD方法來處理所測量的位場數據，由此來分析位場數據所包含的各種成

25、分，以及各成分在原數據中的作用，這種方法可以很好的分析數據中存在的噪聲成分，然后有效的去除這些成分。一、EMD分解一維重力異常數據下列圖為某非線性、非平穩的一維重力數據。圖5 一維重力異常數據它的EMD處理結果如下：圖6 重力數據的EMD分解從圖上可以清楚的看到：重力數據進行EMD分解后，IMF1和IMF2是原始信號中的高頻成分，而IMF3、IMF4和IMF5是低頻成分，在重力異常數據中，高頻的局部代表局部重力異常，而低頻的局部代表區域重力異常。二、EMD分解二維重力異常數據下列圖是某地區非平穩、非線性的重力異常剖面，用EMD方法對其進行分解，我們可以判斷重力異常的類型。圖7 重力異常剖面此重

26、力異常數據的EMD分解如下：圖8 經驗模態分解后IMF1圖象圖9 經驗模態分解后IMF2圖象圖10 經驗模態分解后IMF3圖象比照經EMD處理過的IMF1、IMF2以及IMF3的圖象可以看出以下特征：1IMF1圖象平穩，重力異常值集中，頻率高，表征著重力局部異常。2IMF2圖象的重力異常值明顯比IMF1圖象跨度大，頻率低，是區域異常的特征。3IMF3圖象中，重力異常值變化很大，圖象較IMF2的頻率更低，對研究重力異常的意義不大。 需要說明的是此區域是基于Y軸方向固定，向X軸方向延伸進行數據處理的。而基于X軸固定，數據向Y軸方向延伸的EMD處理結果如下： 圖11 Y方向EMD分解IMF1圖12

27、Y方向EMD分解IMF2由于IMF3和IMF4對原信號的研究意義不大，因此，只分解了2個固有模函數。和以上對重力異常的EMD分解有相同的結論。第二節 EMD分解地震資料地震波形是具有時變特性(或稱非穩態性質)的典型信號,對于這類信號,不僅需要從總體上了解它的頻率成分,而且還需要了解每一時刻信號中所包含的頻率成分。從上文的分析可以看出：經驗模態分解(EMD),能有效地將信號的各種頻率成分以固有模態函數(IMF)形式從時間曲線中別離出來。不同的IMF分量是平穩信號或簡單的非線性信號,具有簡單的非線性特征。其緩變波包特征意味著不同特征尺度波動的波幅隨時間變化,因而也具有時間上的局域化特征,而IMF那

28、么屬于窄帶信號,正好滿足Hilbert變換的要求。對IMF序列進行Hilbert變換,得到包含時間、頻率、振幅的三維離散譜,可提供更加清晰的局部細節時頻特征。一、EMD分解一維地震數據下列圖是某非線性非平穩一維地震數據。圖13 一維地震信號它的EMD處理過程如圖15，分析如下：1原信號共分解出8個固有模函數，其中IMF8缺乏一個周期，為剩余分量。2隨著EMD分解的進行，各IMF分量的頻率逐漸降低，波長逐漸增長，各IMF分量包含了不同的時間特征尺度,可以用不同的分辨率顯示信號特征,并且這種分辨率是根據信號本身的性質自適應的。3隨著模態函數序號的增大,最大振幅根本上呈現逐漸遞減的趨勢,IMF1的最

29、大振幅比其他模態大,說明這一波動所占能量最大。4EMD最先分解出的幾個IMF分量(IMF1IMF4)表達了原始信號中最顯著的信息,是該地震信號的優勢頻段，從圖14可以看出。從C5開始波動的幅值都小于,中心頻率都低于1Hz。依EMD的本性而言,只要有多于1個波的波動存在,EMD就能夠把它提取出來。對于具體的數據資料,這些振幅很小、頻率極低、波長很大的波動可能是事實存在的物理現象,也可能是由于數據采樣率不夠和波動銜接等原因造成的。圖14 IMF1 + + IMF4圖(5)圖15中x2(t)是各個IMF的疊加。和x1(t)相比，根本上沒有失真。這也說明了EMD方法分析信號的可靠性。圖15 一維地震數

30、據的EMD處理過程二、EMD分解地震時間剖面圖16為某地震時間剖面圖，通過EMD方法我們可以更清楚的知道原函數包含的不同成分的信息。圖16 某地震時間剖面圖以下是EMD分解情況。只分解出3個固有模函數圖17 某地震時間剖面的IMF1圖18 某地震時間剖面的IMF2圖19 某地震時間剖面的IMF3從IMF1的圖象中，可以清楚的看到原函數中高頻信息。IMF函數可以進行很好的Hilbert變換，通過對瞬時頻率，瞬時相位，瞬時振幅的分析，可以得到更多的信息。見下章三、EMD分解地震時間切片下列圖為某地震時間切片圖：圖20 地震時間切片圖它的EMD分解如下：圖21 IMF1圖象圖22 IMF2圖象比照地

31、震時間切片圖和經過EMD處理后的IMF1，IMF2的圖象可以得知：原圖象的信息在IMF1中表達的非常明顯，而且輪廓清晰，濾除了一些其他的成分。由于IMF3和IMF4對原函數分析意義不大，不再分析。第三章 基于EMD的Hilbert變換與地震信號屬性信息第一節 基于EMD的Hilbert變換原理進行EMD分解的目的之一就是對信號進行Hilhert變換，進而得到一系列Hilbert譜。基于IMF分量的Hilbert譜可以用二維或三維圖形(振幅、頻率和時間)表示。令第i個IMF為， 3-1在對每一個進行Hilbert變換后，得到一個變換平面內的序列： 3-2其中，P為Cauchy主值。由信號和可以構

32、成一個復數，稱為解析函數。的表達式如下： 3-3這樣，地震屬性信息的瞬時振幅： 3-4瞬時相位： 3-5瞬時頻率： 3-6 3-7因此，原始數據序列可以表示為: 3-8由基于EMD的Hilbert變換根本原理可總結出求Hilbert譜的算法如下：1對每個模函數做Hilbert變換；2原固有模函數和進行Hilbert變換后的模函數根據式3-3組成一個復序列信號；3根據式3-4(3-5)獲得每個時間點的幅值和相位；4按照公式3-6，3-7獲得每一個時間點的瞬時頻率；5繪制瞬時幅值，瞬時相位，瞬時頻率圖譜。第二節 提取地震信號屬性信息地震道在數學意義上是一個解析信號。通過基于EMD的Hilbert變

33、換，我們可以得到地震信號的瞬時振幅譜，瞬時相位譜，瞬時頻率譜。瞬時振幅測定反射強度，它正比例于該時刻地震信號的總能量的均方根值，是識別亮點和暗點的有效工具；瞬時相位是地震剖面上同相軸連續性的量度，相位信息在勾劃有意義的現象像尖滅、斷層、上超、前積反射時很有用；瞬時相位變化的時間速率就是瞬時頻率，瞬時頻率變化很大，它與地層有關，對識別某些地層有幫助，并且會衰減高頻。一、提取地震時間剖面屬性信息下列圖為圖17所示圖形原始地震剖面經EMD處理所得到的固有模態函數IMF1的瞬時振幅，瞬時相位，瞬時頻率圖。圖23 IMF1的Hilbert瞬時振幅圖從IMF1的Hilbert變換的瞬時振幅圖上可以看出：當

34、t=50s,130s時，地震波的反射強度比擬強。圖24 IMF1的Hilbert瞬時相位圖從瞬時相位圖上，可以清楚的顯示地震剖面上同軸連續性。圖25 IMF1的Hilbert瞬時頻率圖從瞬時頻率圖上可以看出從0s到120s之間的頻率與從120s到250s之間的頻率有較大的變化，這與地層有一定關系。二、提取地震時間切片屬性信息下列圖為圖21所示的IMF1經過Hilbert變換以后得出的瞬時振幅圖象：圖26 地震時間切片的IMF1的瞬時振幅圖上可以顯示出時間切片XY平面上各處的幅值，為分析地震波的反射強度提供了有效的信息。結 論本文首先研究了EMD信號分析方法的根本原理和算法，對地球物理資料地震資

35、料，重磁資料進行EMD分解試驗分析，然后研究了基于EMD的Hilbert變換原理及其在提取地震屬性信息的應用，對實際地震時間剖面和時間切片做了試驗分析。本文的方法研究和數據試驗分析說明：1、經EMD分解變換得到的IMF序列是直接從原始時序數據中別離出來的,事先無需確定分解階次,不存在機械分解。因此IMF序列能更好反映原始數據固有的物理特性,其分解是客觀的,內在的和自適應的。每階IMF序列都代表了某種特定意義的頻帶信息,給實際應用與解釋工作帶來了方便；2、原始信號經過EMD分解獲得的IMF序列具有穩態性,對IMF進行Hilbert變換，就可以得到單個固有模態函數的瞬時振幅圖，瞬時頻率圖，瞬時相位

36、圖，這些圖可以清楚的顯示信號的時頻特征；3、EMD分析方法用于分解地球物理資料和作時頻分析是有效的。下一步的工作是深入研究地球物理信號EMD分解結果的地球物理意義。致 謝本次畢業設計得到了孟小紅教授的悉心指導以及深切的關心，孟老師在百忙之中能夠親自指導開題報告和中期報告以及最后的審閱工作，在此表示深深的感謝和崇高的敬意。本次畢業設計從始至終都得到了郭良輝老師的具體指導，包括對論文的修改工作都做到了無微不至。郭老師知識出眾，對待同學和藹可親，對待知識嚴肅認真，在很多方面都給我留下了深刻的印象。在此真誠的向郭老師道一聲：您辛苦了，謝謝。在本科期間，得到了諸多任課老師的培養，幫助和關心，這里表示由衷

37、的感謝。與同學們四年的相處也讓我學到了很多知識，得到了很多快樂，在此一并表示真心的感謝。計算機專業畢業設計開發環境：ASP.NET，VB，VB.NET VF，java等，數據庫：SQL。包括:開題報告、程序、論文、辯論PPT，所有程序都是通過辯論的優秀作品，質量保證。也可代做。我是哈爾濱工業大學計算機專業畢業的學生我賣的畢業設計都是新做出來的而且是學生親手做的符合學生要求如果你在別的店賣來的都是很專業的人員做的一看就不是學生自己親手做出來的，而且其他店不提供售后我們提供售后效勞及技術支持和辯論技巧.Q Q：982465840本店購設計的優點: 1價格合理廉價2提供技術支持3售后效勞好4成交速度

38、快當時就可以完成調試功能5東西齊全(開題論文代碼程序辯論PPT售后效勞)6作品都是獲得優秀的產品(保證質量)參考文獻1 程佩青. 數字信號處理. 北京:清華大學出版社,1994.2 Huang, N.E., Shen, Z., Long, S.R., et. al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis. Royal Society London, 1998: 903-995.3 熊學軍, 郭炳火. EMD方法和

39、Hilbert 譜分析法的應用與探討. 黃渤海海洋, 2002, 20(2). 12-21.4 武安緒, 吳培稚, 等. Hilbert-Huang 變換與地震信號的時頻分析. 中國地震, 2005, 21(2): 207-215.5 Hassan H. Hassan, GEDCO. Empirical Mode Decomposition (EMD) of potential field data. SEG 72th Annual International Meeting Expanded Abstracts, 2005: 704-706.6 胡勁松. 面向旋轉械故障診斷的經驗模態分解時頻

40、分析方法及實驗研究. 博士學位論文. 杭州： 浙江大學, 2003.7 朱金龍. 經驗模態分解方法及在癲癇腦電波中的應用研究. 碩士學位論文. 南京: 南京郵電大學. 2006.8 秦勃, 時鵬. 經驗模式多尺度圖象分解. 工程圖學學報, 2002, 4: 73-78.9 張志涌等編著. 精通Matlab6.5版. 北京: 北京航空航天大學出版社, 1998.10 管志寧. 地磁場與磁力勘探. 北京:地質出版社,2005.11 曾華霖. 重力場與重力勘探. 北京:地質出版社,2005.12 羅孝寬, 郭紹雍. 應用地球物理教程重力、磁法. 北京: 地質出版社, 1991.附錄：經驗模態分解MA

41、TLAB程序經驗模態分解程序function imf,ort,nbits = emd(varargin);x,t,sd,sd2,tol,display_sifting,sdt,sd2t,ner,nzr,lx,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask = init(varargin:);if display_sifting figureend% maximum number of iterations% MAXITERATIONS=2000;%main loop : requires at least 3 e

42、xtrema to proceedwhile stop_EMD(r) & (k MAXMODES+1 | MAXMODES = 0) & any(mask)% current modem = r;% mode at previous iterationmp = m;if FIXEstop_sift,moyenne = stop_sifting_fixe(t,m,INTERP);elseif FIXE_Hstop_count = 0;stop_sift,moyenne,stop_count = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H);s

43、top_count = 0;elsestop_sift,moyenne = stop_sifting(m,t,sd,sd2,tol,INTERP);end if (max(m) - min(m) (1e-10)*(max(x) - min(x) if stop_sift warning(forced stop of EMD : too small amplitude) else disp(forced stop of EMD : too small amplitude) end break end % sifting loop while stop_sift & nbitMAXITERATIO

44、NS/5 & mod(nbit,floor(MAXITERATIONS/10)=0 & FIXE & nbit 100) disp(mode ,int2str(k), iteration ,int2str(nbit)if exist(s) disp(stop parameter mean value : ,num2str(s)end im,iM = extr(m); disp(int2str(sum(m(im) 0), minima 0; ,int2str(sum(m(iM) 0), maxima 100) if exist(s) warning(forced stop of sifting

45、: too many iterations. mode ,int2str(k),. stop parameter mean value : ,num2str(s) else warning(forced stop of sifting : too many iterations. mode ,int2str(k),.) end end end % sifting(篩選) loop imf(k,:) = m; if display_sifting disp(mode ,int2str(k), stored) end nbits(k) = nbit; k = k+1; r = r - m; nbi

46、t=0;end %main loopif sum(r.2) & any(mask)imf(k,:) = r;endort = io(x,imf);if display_sifting closeend%function stop = stop_EMD(r)%停止EMDindmin,indmax,indzer = extr(r);ner = length(indmin) + length(indmax);stop = ner sd) tol | any(sx sd2) | (abs(nzm-nem)1) & (nem 2);catchdisp(lasterr)stop = 1;envmoy =

47、zeros(1,length(m);% disp(catch : ,lasterr)s = NaN;end%function stop,moyenne= stop_sifting_fixe(t,m,INTERP)tryenvmin,envmax,moyenne = envelope(t,m,INTERP);stop = 0;% disp(try ok)catchmoyenne = zeros(1,length(m);stop = 1;% disp(catch : ,lasterr)%lasterr:matlab發出的最新錯誤信息end%function stop,moyenne,stop_co

48、unt= stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H)tryenvmin,envmax,moyenne,indmin,indmax,indzer = envelope(t,m,INTERP); nem = length(indmin) + length(indmax); nzm = length(indzer);if (abs(nzm-nem)1)stop = 0;stop_count = 0;elsestop_count = stop_count+1;stop = (stop_count = FIXE_H);end% disp(try o

49、k)catchmoyenne = zeros(1,length(m);stop = 1;% disp(catch : ,lasterr)end%function display_emd(t,m,mp,r,envmin,envmax,envmoy,s,sb,sx,sdt,sd2t,nbit,k,display_sifting,stop_sift)subplot(4,1,1) plot(t,mp);hold on; plot(t,envmax,-k);plot(t,envmin,-k);plot(t,envmoy,r); title(IMF ,int2str(k),; iteration ,int

50、2str(nbit), before sifting); set(gca,XTick,) hold off subplot(4,1,2) plot(t,sx) hold on plot(t,sdt,-r) plot(t,sd2t,:k) title(stop parameter) set(gca,XTick,) hold off subplot(4,1,3) plot(t,m) title(IMF ,int2str(k),; iteration ,int2str(nbit), after sifting); set(gca,XTick,) subplot(4,1,4); plot(t,r-m)

51、 title(residue); disp(stop parameter mean value : ,num2str(sb), before sifting and ,num2str(s), after)if stop_siftdisp(last iteration for this mode)end if display_sifting = 2 pause(0.01) else pause end%function display_emd_fixe(t,m,mp,r,envmin,envmax,envmoy,nbit,k,display_sifting)subplot(3,1,1) plot

52、(t,mp);hold on; plot(t,envmax,-k);plot(t,envmin,-k);plot(t,envmoy,r); title(IMF ,int2str(k),; iteration ,int2str(nbit), before sifting); set(gca,XTick,) hold off subplot(3,1,2) plot(t,m) title(IMF ,int2str(k),; iteration ,int2str(nbit), after sifting); set(gca,XTick,) subplot(3,1,3); plot(t,r-m) tit

53、le(residue); if display_sifting = 2 pause(0.01) else pause end %function envmin, envmax,envmoy,indmin,indmax,indzer = envelope(t,x,INTERP)%computes envelopes and mean with various interpolationsNBSYM = 2;DEF_INTERP = spline;if nargin 1error(x and t must be vectors)ends = size(x);if s(1) 1x = x;ends

54、= size(t);if s(1) 1t = t;endif length(t) = length(x)error(x and t must have the same length)endlx = length(x);indmin,indmax,indzer = extr(x,t); %boundary conditions for interpolationtmin,tmax,xmin,xmax = boundary_conditions(indmin,indmax,t,x,NBSYM);% definition of envelopes from interpolationenvmax

55、= interp1(tmax,xmax,t,INTERP);envmin = interp1(tmin,xmin,t,INTERP);if nargout 2 envmoy = (envmax + envmin)/2;end%function tmin,tmax,xmin,xmax = boundary_conditions(indmin,indmax,t,x,nbsym)% computes the boundary conditions for interpolation (mainly mirror symmetry)lx = length(x);if (length(indmin) +

56、 length(indmax) 3)error(not enough extrema)endif indmax(1) x(indmin(1)lmax = fliplr(indmax(2:min(end,nbsym+1);lmin = fliplr(indmin(1:min(end,nbsym);lsym = indmax(1);elselmax = fliplr(indmax(1:min(end,nbsym);lmin = fliplr(indmin(1:min(end,nbsym-1),1;lsym = 1;endelseif x(1) x(indmax(1)lmax = fliplr(in

57、dmax(1:min(end,nbsym);lmin = fliplr(indmin(2:min(end,nbsym+1);lsym = indmin(1);elselmax = fliplr(indmax(1:min(end,nbsym-1),1;lmin = fliplr(indmin(1:min(end,nbsym);lsym = 1;endend if indmax(end) indmin(end)if x(end) x(indmin(end)rmax = fliplr(indmax(max(end-nbsym,1):end-1);rmin = fliplr(indmin(max(en

58、d-nbsym+1,1):end);rsym = indmax(end);elsermax = fliplr(indmax(max(end-nbsym+1,1):end);rmin = lx,fliplr(indmin(max(end-nbsym+2,1):end);rsym = lx;endend tlmin = 2*t(lsym)-t(lmin);tlmax = 2*t(lsym)-t(lmax);trmin = 2*t(rsym)-t(rmin);trmax = 2*t(rsym)-t(rmax); % in case symmetrized parts do not extend en

59、oughif tlmin(1) t(1) | tlmax(1) t(1)if lsym = indmax(1)lmax = fliplr(indmax(1:min(end,nbsym);elselmin = fliplr(indmin(1:min(end,nbsym);endif lsym = 1error(bug)endlsym = 1;tlmin = 2*t(lsym)-t(lmin);tlmax = 2*t(lsym)-t(lmax);end if trmin(end) t(lx) | trmax(end) 2x1=x(1:m-1);x2=x(2:m);indzer = find(x1.

60、*x20);if any(x = 0) iz = find( x=0 ); indz = ; if any(diff(iz)=1) zer = x = 0; dz = diff(0 zer 0); debz = find(dz = 1); finz = find(dz = -1)-1; indz = round(debz+finz)/2); else indz = iz; end indzer = sort(indzer indz);endend d = diff(x);n = length(d);d1 = d(1:n-1);d2 = d(2:n);indmin = find(d1.*d20