趙應權， 沈忠民， 周愛芬

(1. 成都理工大學 a.“油氣藏地質及開發(fā)工程”國家重點實驗室，b.能源學院，成都 610059;2．中國石油工程建設公司，北京 100120)

?

Mallat小波濾波器系數(shù)的計算與地質應用

趙應權1a，b， 沈忠民1a， 周愛芬2

(1. 成都理工大學 a.“油氣藏地質及開發(fā)工程”國家重點實驗室，b.能源學院，成都 610059;2．中國石油工程建設公司，北京 100120)

針對第一類、第二類Mallat系列小波的濾波器公式為頻率域的表達式，它們無法直接求得對應系列小波的濾波器系數(shù)，從而造成了Mallat系列小波在時間域內多尺度計算的困難。這里基于傅里葉分析原理，借用MATLAB工作平臺編制了由頻率域到時間域的Mallat系列小波濾波器組的求解代碼。該段代碼解決了Mallat系列小波在時間域內的濾波器系數(shù)求解。經實際數(shù)字計算表明，采用該方法求得的濾波器組系數(shù)進行多尺度分析，能夠實現(xiàn)小波的完全分解與重構，這對選用Mallat系列小波中的具體小波進行異常檢測分析，具有十分重要的意義。

Mallat小波； 異常檢測； FFT； 濾波器系數(shù)

0 引言

小波在奇異性信號的檢測上具有獨特的優(yōu)勢，在圖像邊緣檢測、網絡異常檢測、地震異常與語音異常識別等領域有著廣泛地應用[1-5]。然而用不同的小波處理相同的信號數(shù)據(jù)序列有著不同的處理結果，為了取得更好的處理結果， Mallat[6]針對圖像邊緣提取發(fā)明了一類B樣條小波，即Mallat小波。前人利用Mallat小波做過重力異常多尺度反演研究[7]，也用該小波提取過深層次物化探的異常信息[8]，Zhang et.al[9-11]應用該小波在油氣化探異常信號的分離中也取得了良好效果。

Mallat小波分奇數(shù)階B樣條小波和偶數(shù)階B樣條小波兩類。Mallat僅給出了任意階第一類小波，即小波函數(shù)關于原點反對稱的小波函數(shù)ψ、尺度函數(shù)φ、濾波器組H、G、K、L的頻率響應表式。劉曙光[12]補充了任意階第二類B樣條小波，即小波函數(shù)關于原點對稱的小波的推導公式，并給出了該類小波濾波器組H、G、K、L的解析表達式。

由于Mallat小波的濾波器系數(shù)在實際應用中受待處理數(shù)據(jù)影響，需要選用Mallat系列小波中不同階次的小波。因此，不同階次的Mallat小波的濾波器系數(shù)在實際計算中顯得尤為重要。但是在前人研究中，僅有Mallat小波濾波器系數(shù)的頻率域解析公式，工程技術人員很難利用這些公式直接求得不同階次小波的濾波器系數(shù)，這勢必影響到Mallat小波的多尺度計算，從而大大降低了Mallat小波的實用價值。為此作者用Matlab語言求解了Mallat小波的濾波器系數(shù)。

1 濾波器系數(shù)公式

濾波器是實現(xiàn)小波多尺度分析計算的重要參數(shù)。Mallat小波以樣條函數(shù)為平滑函數(shù)θ，以樣條函數(shù)的一階導函數(shù)ψ′(x)或二階導函數(shù)ψ″(x)為小波函數(shù)。Mallat小波的平滑函數(shù)階數(shù)的奇偶性與小波函數(shù)的對稱性不同，濾波器系數(shù)的求解公式也有所不同。具體的求解公式為式(1)與式(2)。

(1)

(2)

在式(1)中，只要平滑函數(shù)θ的階數(shù)相同，對應的第一類小波與第二類小波在小波的分解與重構中有著相同的低通濾波器系數(shù)。在式(2)中，平滑函數(shù)θ的階數(shù)并不對高通濾波器產生作用，僅受小波類型的影響，對于第一類小波，G(ω)為奇函數(shù)，具有關于原點反對稱的特點，對于第二類小波，G(ω)為偶函數(shù)，具有關于原點對稱的特點。這兩類小波的對稱性的不同，正好對應于圖像邊緣檢測的極值檢測與零點檢測兩種檢測方式。

由于Mallat小波分析中的另外兩個濾波器K、L與H、G有著定量的解析關系，它們之間的相關關系為：

(3)

(4)

式(3)、式(4)中的H(ω)、G(ω)對應于式(1)、式(2)中的H(ω)、G(ω)，將參數(shù)H(ω)、G(ω)代入式(3)、式(4)后經傅立葉逆變換可以解出時間域內的K、L濾波器系數(shù)。

2 濾波器系數(shù)求解

根據(jù)前文的濾波器函數(shù)在頻率域的解析表達式，采用快速傅里葉變換(FFT)及逆變換(IFFT)的方法求取Mallat小波的濾波器系數(shù)。具體求解的步驟(圖1)：

1)針對單個頻譜周期，根據(jù)Nyquist采樣定理初步確定采樣點數(shù)N1。

2)針對傅里葉分析要求，修正采樣點數(shù)為N2，且N1≤N2,N2=2N,N∈Z+。

3)根據(jù)修正的采樣點數(shù)N2，求得采樣的平均間隔dw。

4)對一個完整的頻率周期2π按照dw的采樣間隔平均采樣，并計算出各個相位上對應的離散信號序列x[i]。

5)對步驟4)中的離散信號序列x[i]進行傅里葉逆變換和中心移位，得到時間域內的信號序列y[i] 。

6)取步驟5)中的離散信號序列y[i]的模得到實信號序列z[i]，截除實信號序列z[i]兩側的“0”值，剩余的數(shù)據(jù)即為Mallat小波在時間域內的濾波器系數(shù)。

根據(jù)上述操作步驟，采用matlab語言實現(xiàn)了這一求解過程。以平滑函數(shù)θ階數(shù)為奇數(shù)的第一類Mallat小波為例，設定平滑函數(shù)的階數(shù)m=3時，濾波器H、G、K、L的matlab數(shù)字求解代碼如表1所示。

其他濾波器系數(shù)可參考公式(1)-公式(4)用類似的方法求解，運行上述代碼本文求得平滑函數(shù)是3階、4階Mallat小波的濾波器系數(shù)如表2所示。

3 應用實例

3.1 地表油氣化探干擾識別

東營凹陷金東-柳橋地區(qū)地表油氣化探數(shù)據(jù)中干擾嚴重，為從這些數(shù)據(jù)中提取有效的異常信息，需要排除干擾。前人基于異常與干擾的頻率差異，發(fā)明了小波排除干擾的方法[9-11]。針對金東-柳橋地區(qū)的地表油氣化探數(shù)據(jù)，作者采用平滑函數(shù)階次為3的Mallat反對稱小波，進行了干擾排除與異常信息的識別。圖2為該區(qū)酸解烴甲烷在東西測線上的干擾與異常分析圖，其中低頻序列對應著干擾信息，高頻序列對應著異常信息，由圖2可以看出，無論是干擾，還是異常在空間位置上都一一對應，并且具有繼承性，這為后續(xù)的干擾排除，準確地提取異常信息提供了可能。

圖1 Mallat小波濾波器系數(shù)數(shù)字求解過程圖

表1 第一類Mallat小波的濾波器系數(shù)求解代碼

3.2 層序地層劃分

元壩地區(qū)是中國石化繼普光氣田后在川東北地區(qū)臺地邊緣超深層領域內，又一個取得重大突破的區(qū)塊。前人的大量研究認為，川東北地區(qū)長興組儲層分布明顯受Ⅲ級層序體系域及沉積相控制，儲層主要分布于SQ1、SQ2及SQ3的高位體系域中[13]。因而層序的快速準確劃分，為長興組—飛仙關組儲層的時空展布提供預測依據(jù)，也可促進該區(qū)天然氣勘探的快速發(fā)展。

地層縱向的周期性為使用小波劃分層序提供了理論依據(jù)，前人在實踐中發(fā)現(xiàn)不同小波在層序劃分中存在較大差異[14]，即小波分析手段劃分出的層序界面位置圍繞實際層序界面位置的上下波動，因此采用小波方法劃分層序往往需要大量的人工調整。考查這些用于劃分層序的小波發(fā)現(xiàn)，前人使用的小波大都為Daubechies小波、Symlets小波和DMeyer小波，這些小波在對稱性或緊支性上較差，容易造成分析數(shù)據(jù)序列中的異常值的“漂移”，進而導致層序界面點不準。而Mallat小波作為邊緣檢測的最佳小波，正好客服了上述不足。以元壩地區(qū)的122井二疊統(tǒng)的長興組地層為例，結合地層標志人工劃分該段地層為2個三級層序，7個四級層序，16個五級層序(圖3)。采用小波劃分該段地層時發(fā)現(xiàn)，DB4小波劃分出的層序界面位置較人工解釋界面位置有較大偏差，相鄰的四級、三級層序在層序界面點上也并不具有很好的繼承性；而用平滑函數(shù)階次為3的Mallat反對稱處理的結果跟人工解釋的層序界面完全一致，這樣就大大減少了人工調整層序界面位置的工作量。

表2 平滑函數(shù)階數(shù)為3、4的Mallat小波濾波器系數(shù)

注：m.平滑函數(shù)的階數(shù)，T.小波類型，H.低通濾波器，G.高通濾波器，K.重構高通濾波器，L.二維小波重構濾波器

圖2 金東-柳橋地區(qū)酸解烴甲烷東西測線上的Mallat小波分析圖

圖3 川東北元壩122井地層層序劃分圖

4 結論

基于傅里葉分析原理，利用Matlab平臺求解了Mallat系列小波在時間域內的濾波器系數(shù)。傅里葉分析是聯(lián)接時間域與空間域計算的紐帶，Matlab軟件的數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計算簡潔便利、人機交互能力強的特點，為Mallat小波系數(shù)的求解提供了條件。在Nyquist采樣定理的約束下，對周期范圍內的數(shù)據(jù)實行采樣離散化，然后利用傅里葉分析手段，在Matlab工作平臺編制的由頻率域到時間域的Mallat系列小波濾波器組的求解代碼，這組代碼解決了mallat系列小波在時間域內的濾波器系數(shù)求解。

[1] 韓慧妍，韓燮．基于方向小波變換的邊緣檢測算法[J]. 微電子學與計算機，2012, 29(7): 55-57. HAN H Y, HAN X. Edge detection algorithm based on directional wavelet transform[J]. Microelectronics & computer. 2012, 29(7): 55-57.(In Chinese)

[2]馮興杰，焦文歡. 基于動態(tài)閾值的網絡異常檢測[J]. 計算機工程與設計，2012, 33(6): 2182-2186. FENG X J, JIAO W H. Network abnormality detection based on dynamic threshold[J]. Computer engineering and design, 2012, 33(6): 2182-2186.(In Chinese)

[3]解滔，杜學彬，劉君，等．汶川Ms8.0、海地Mw7.0地震電磁信號小波能譜分析[J]．地震學報，2013, 35(1):61-71. XIE T, DU X B, LIU J, et al. Wavelet power spectrum analysis of the electromagnetic signals of Wenchuan Ms8.0 and Haiti Mw7.0 earthquakes[J]. Acta seismologicaca sinica. 2013, 35(1):61-71.(In Chinese)

[4]李衛(wèi)，宋弘，姜天華．基于小波神經網絡的嵌入式語音識別系統(tǒng)[J]. 通信技術, 2010, 43(6): 213-218. LI W, SONG H, JIANG T H. Embedded speech recognition system based on wavelet neural network[J]. Communications Technology. 2010, 43(6): 213-218.(In Chinese)

[5]敬少群，劉春平, 王佳衛(wèi)．小波變換在深井水位異常識別中的應用[J]. 地震研究, 2012, 35(2): 171-176. JING S Q，LIU CH P，WANG J W. Application of wavelet multi-scaling decomposition in identifying water level anomaly of deep well[J]. Journal of seismological research. 2012, 35(2): 171-176.(In Chinese)

[6]STEPHANE MALLAT, SIFEN ZHONG. Characterization of signals from multiscale edges[J]. IEEE Transactions on pattern analysis and machine intelligence, 1992, 14(7): 710-732.

[7]LIUPING ZHANG, GUOPING BAI,KEBIN ZHAO. Restudy of acid-extractable hydrocarbon data from surface geochemical survey in the Yimeng Uplift of the Ordos Basin, China: Improvement of geochemical prospecting for hydrocarbons[J]. Marine and Petroleum Geology, 2006(23): 529-542.

[8]LIUPING ZHANG, GUOPING BAI. A wavelet-analysis-based new approach for interference elimination in geochemical hydrocarbon exploration[J]. Mathematical Geology, 2003, 35(8): 939-952.

[9]ZHANG LIUPING, RUAN TIANJIAN. Application of wavelet analysis to interference elimination for geochemical hydrocarbon exploration[J]. Journal of China University of Geosciences, 2000, 11(1): 89-91.

[10]侯遵澤，楊文采．重力異常多尺度反演研究及應用[C]．中國地球物理學會年會，2011. HOU Z, YANG W C. Study on multi-scale inversion of gravity anomalies[C]. The Chinese Geographics annual meeting, 2011.(In Chinese)

[11]陳建國，夏慶霖．利用小波分析提取深層次物化探異常信息[J]．地球科學-中國地質大學學報，1999， 24(5)：509-512. CHEN J G, XIA Q L. Wavelet-based extraction of geophysical and geochemincal anomaly information[J]. Earth Science—Journal of China University of Geosciences, 1999, 24(5): 509-512.(In Chinese)

[12]劉曙光．基于邊緣檢測算子的兩類B樣條小波[J]．紡織高?；A科學學報，2001，14(1)：66-71. LIU SH G. Two types of the B-spline wavelets based on edge detector[J]. Basic sciences journal of textile universities. 2001, 14(1): 66-71.(In Chinese)

[13]郭彤樓．川東北元壩地區(qū)長興組—飛仙關組臺地邊緣層序地層及其對儲層的控制[J]．石油學報，2011，32(3)：387-394. GUO T L. Sequence strata of the platform edge in the Changxing and Feixianguan formations in the Yuanba area,northeastern Sichuan basin and their control on reservoirs[J]. Acta petrolei sinica, 2011, 32(3): 387-394.(In Chinese)

[14]朱劍兵，紀友亮，趙培坤，等．小波變換在層序地層單元自動劃分中的應用[J]．石油勘探與開發(fā)，2005，32(1)：84-86. ZHU J B, JI Y L, ZHAO P K, et al. Application of wavelet transform in auto-identify units of stratigraphy sequence[J]. Petroleum Exploration and Development, 2005, 32(1):84-86.(In Chinese)

Geological application and calculation of Mallat wavelet filter coefficients

ZHAO Ying-quan1a,b, SHEN Zhong-min1a, ZHOU Ai-fen2

(1.Chengdu University of Technology a. State Key Laboratory of Oil & Gas Reservoir Geology and Exploration, b. college of Sedimentary，Chengdu 610059, China；2.China Petroleum Engineering & Construction Corp., Beijing 100120,China)

Formulae for the 1st and 2nd type of Mallat series wavelet filters are expressed in frequency domain, the wavelet filter coefficients cannot obtained by the corresponding formula directly. Thus, it is difficult to calculate in multi-scale analyzing for the Mallat series wavelet in time domain. In this paper, based on the principle of Fourier analysis, some solving code is complied for the Mallat series wavelet filter's from the frequency domain to the time domain by the MATLAB language. This code solves the filter coefficients for Mallat series of wavelet in time domain. The actual number calculation shows that decomposition and reconstruction of wavelet calculate in multi-scales can be achieved completely by the filter coefficients from the method above-mentioned, and it has great significance for anomaly detection to select a specific Mallat wavelet.

Mallat wavelet; anomaly detection; FFT; filter coefficient

2013-12-23 改回日期：2014-09-15

國家科技重大專項項目(2008ZX05008-004-20);“油氣藏地質及開發(fā)工程”國家重點實驗室開放基金項目( PLC201102)

趙應權(1977-)，男，博士，主要從事油氣勘探研究工作，E-mail: hbszzyq@126.com。

1001-1749(2015)04-0532-06

P 628

A

10.3969/j.issn.1001-1749.2015.04.19