徐彥凱,雙 凱
(中國石油大學地球物理與信息工程學院,北京 102249)
在油氣井下電磁脈沖數據傳輸中,有效信號具有暫態(tài)瞬變的特征。由于傳輸距離遠,接收到的微弱有用信號往往被環(huán)境噪聲所淹沒。因此,有效檢測瞬變弱信號是電磁脈沖數據傳輸的關鍵技術。小波分析是近些年發(fā)展起來的一種新的時頻分析方法,小波變換能使信號能量集中在一些大的小波系數中;而噪聲能量卻分布于整個小波域內,對應的小波系數小。因此,小波閾值法可以較好地用于抑制噪聲并檢測信號??紤]到加噪信號對應矩陣的奇異值分解一定程度上可以反映矩陣的主要特征,實現噪聲去除和信號提取,筆者將小波及奇異值分解相結合,實現微弱瞬變信號的檢測。
設ψ(t)∈L2(R),若ψ(t)為母小波,其傅里葉變換Ψ(ω)滿足條件:
將ψ(t)進行伸縮和平移,得到小波基函數
其中a為尺度(伸縮)因子,b為平移因子。信號f(t)的連續(xù)小波變換定義[1]為
對 a、b 進行二進離散化,即 a=2j,b=2jk,j,k∈Z,得到的二進小波變換為
其逆變換為
Mallat算法實現了二進小波的快速分解與重構[2],分解公式為
重構公式為
式中,cj,k為尺度系數;dj,k為小波系數;h、g、~h、~g 分別為低通和高通分解與重建濾波器;j為分解層數。
由于小波變換能使信號的能量集中在一些大的小波系數中,而噪聲能量卻平均分布于整個小波域內,對應的小波系數小。應用信號與噪聲小波系數的特點對信號降噪,降噪步驟[3]如下:
(1)確定小波基函數和小波分解最大層數,應用式(5)對信號進行J層小波分解得到小波系數和尺度系數。
(2)確定每層小波系數的閾值,根據每層閾值按閾值函數形式對小波系數進行處理,把保留下來的作為降噪后信號的小波系數。
(3)利用式(6)對信號進行重構,并對降噪后信號檢測。
上述步驟中,分解尺度、閾值函數和小波基函數的合理選擇對信號檢測至關重要。
2j尺度下小波偽頻率f2j與小波中心頻率fc的關系[3]:
式中,Δt為采樣周期。
若信號的中心頻率接近小波的偽頻率,則會使此頻率段的小波系數達到最大,從而有效地檢測信號。研究的瞬變信號的頻率集中在(0.3~1.1)×105Hz,db2小波的中心頻率為0.6667 Hz,采樣周期為1×10-6s,由式(7)可知,j=3和4的db2小波偽頻率在信號的頻率范圍,所以確定小波分解最大層數為4。
常用的方法有硬閾值函數和軟閾值函數,采用軟閾值處理會對所有的小波系數都進行抑制,可使降噪后信號比較平滑,但對于提取具有突變特性的瞬變信號而言,對降噪后信號的平滑性沒有要求,相反采用硬閾值方法能夠更好地保留原信號的一些瞬變特征以便信號檢測,所以采用硬閾值處理方法更為合適。
選擇小波基函數時考慮正交小波變換后的數據不會出現冗余,并保證信號精確重構。因此選擇幾種有代表性的常用正交小波基:haar小波、coif1小波、sym2小波和db2小波,在此基礎上進一步通過仿真試驗對比選優(yōu)。
四種小波基下的重建信號見圖1。
圖1 四種小波基下的重建信號Fig.1 Reconstructed waveforms by four wavelet functiones
電磁脈沖傳輸的數字信號用圖1(a)中的瞬變信號表示,該信號是一個持續(xù)時間短的暫態(tài)過程。通過電路分析,得到該信號的數學模型為
式中,A為信號幅度;ω為振蕩頻率;α為衰減因子,決定信號暫態(tài)過程的持續(xù)時間,其持續(xù)時間由電路參數決定;t≤t0為充電時間;t>t0為放電時間。α=15,ω =4.7 × 105rad/s。
井下電磁脈沖信號傳輸中,噪聲源主要有兩大類:井下鉆井環(huán)境產生的噪聲和電子器件產生的電子噪聲。鉆井噪聲頻率范圍為1~4 kHz。電子噪聲在頻域和時域上分布一致,是一種高斯白噪聲[4]。由于鉆井噪聲可通過濾波的方法去噪,在此只考慮無法直接濾除的白噪聲。
針對小波基的選取進行試驗,試驗中采用如圖1(b)所示的-2.5 dB加噪信號,統(tǒng)一的4尺度分解,得到圖1所示的4種小波基下的重構信號。比較圖1中的(c)、(d)、(e)、(f),可以直觀地看出,針對本文研究的瞬變信號,db2小波的降噪效果最佳,sym2次之,coif1小波最差,因此選擇db2小波基。
由圖1(b)信號作4尺度db2小波分解的各尺度系數如圖2所示。從圖2看出,無法從j=1,2的小波系數中區(qū)分信號和噪聲,因此重構時將j=1,2的小波系數全置零,保留 j=3,4中大于閾值 λ=σ的小波系數。
圖2 db2小波分解各尺度系數Fig.2 Scale coefficients with db2 wavelet
圖1(c)表明,在高信噪比的情況下,小波閾值法可以實現信號的檢測,但從圖3所示的 -8.5 dB加噪信號的db2小波重建信號可以看出:在(0.1~0.2)ms和(0.7~0.8)ms處沒有信號,但在相應位置重構信號的幅度較大。
因此,為了更好地檢測強噪聲中的弱信號,需要將j=3,4信號和噪聲的小波系數進一步分離。
圖3 -8.5 dB加噪信號的db2小波重建信號Fig.3 db2 WT reconstructed waveform with-8.5 dB SNR
設矩陣Hm×n(m≤n)秩為r,則存在兩個正交矩陣U、V和對角矩陣D,使得下式成立[8]:
式中,λi為矩陣H的奇異值;ui和vi是λi的特征向量為特征項。奇異值分解是一種代數特征提取方法,能找到矩陣各向量之間的內在本質屬性。因此加噪信號序列的奇異值分解可以提取序列主要特征,實現信號降噪和檢測。其步驟如下:
(1)Hankle矩陣生成。設加噪信號序列為X=[x0,x1,…,xN-1],構造 Hankle 矩陣 H[5-7]
矩陣H的行數L、列數J及序列X的長度N之間應滿足
(2)矩陣奇異值分解。對矩陣H進行奇異值分解,λ1,λ2,…,λr為矩陣 H 的按降序排列的非零奇異值,即 λ1≥ λ2≥ … ≥ λr。
由式(12)重建信號并檢測有用信號存在與否。
奇異值分解中,窗長L及保留特征值個數TR的選擇影響降噪和檢測效果。
本文首先選取TR=5,窗長L分別選取15、22、30、40、60,碼元中有用信號的采樣點數是 30,加噪信號的信噪比分別為-2.5、-5.5、-8.5、-9.5、-10.5、-11.5 dB,仿真得到如圖4所示的SVD降噪峰值信噪比與窗長L的關系。
圖4 窗長對SVD降噪峰值信噪比的影響Fig.4 Effect of window length on SVD-denoised PSNR
從圖4看出:當原信噪比大于-8.5 dB,窗長L選取30效果最好;當原信噪比小于-8.5 dB,窗長L選取22效果最好;窗長L選取過大(L=60)或過小(L=15)降噪效果明顯差。因此得出:①若原信噪比較大,選擇窗長L等于碼元中有用信號長度時效果最佳;②若原信噪比較小,選擇窗長L比碼元中有用信號的長度稍小些降噪和檢測效果最好,這是由于信號末端采樣值接近零,當噪聲較強時,奇異值分解中這部分值的特征完全被噪聲淹沒造成的。
然后選取窗長L=22,保留特征值個數TR分別選取15、10、5、3、2,加噪信號的信噪比分別為 -2.5、-5.5、-8.5、-9.5、-10.5、-11.5 dB,得到如圖5所示的SVD降噪峰值信噪比與保留特征值個數TR的關系。
圖5 保留特征值個數TR對SVD降噪峰值信噪比的影響Fig.5 Effect of main singular values1 number on SVD-denoised PSNR
從圖5看出:①降噪后峰值信噪比隨著原信噪比的減小而減小;②保留特征值個數TR=2峰值信噪比最高,檢測效果也最佳;③原信噪比越小,TR的選擇對降噪后峰值信噪比的影響越大,這是由于奇異值分解中大部分特征值已經體現不出信號的特性。
在研究奇異值分解參數選擇的基礎上,采用奇異值分解方法對-8.5 dB加噪信號檢測,其中選擇參數L=22,TR=2,仿真結果如圖6所示。
圖6 -8.5 dB加噪信號的SVD重建信號Fig.6 Reconstructed waveform by SVD method with-8.5 dB SNR
比較圖3(c)和圖6(c)可知:從檢測的角度看,奇異值分解比小波降噪后判決出錯的概率小;但從輸出信噪比看,奇異值分解比小波降噪差。
通過上述分析可知,若將小波和奇異值相結合[7-8]會使信號和噪聲的小波系數更好地分離,檢測效果更佳。具體方法是:首先將信號作小波分解,再對小波分解系數作奇異值分解,最后通過閾值法保留信號小波系數并重建降噪信號,利用重建信號進行信號檢測。
圖7(a)和圖7(b)分別為信號和-8.5 dB加噪信號,采用小波閾值、奇異值分解及小波奇異值三種方法對圖7(b)加噪信號降噪,降噪結果分別如圖7(c)~(e)??梢钥闯?僅采用奇異值分解的圖7(d)降噪效果最差;采用小波降噪的圖7(c)雖然效果較好,但在沒有信號的(0.1~0.4)ms、和(0.5~0.6)ms處的重構信號的幅度較大;采用小波奇異值方法的圖7(e)則可以準確檢測信號的存在。因此與小波閾值法和奇異值法相比,小波奇異值法降噪和檢測效果最佳。
圖7 -8.5 dB加噪信號的三種方法重建信號Fig.7 Reconstructed waveforms by three methods with-8.5 dB SNR
針對油氣井下電磁脈沖數據傳輸中瞬變弱信號的檢測,研究了小波閾值法和奇異值分解及參數選擇。在此基礎上提出了小波與奇異值分解相結合降噪檢測信號的方法。仿真結果表明,該方法能更好地區(qū)分信號和噪聲,獲得了更好的降噪和檢測結果。
[1] IMTIAZ H,FATTAH S A.A wavelet-based dominant feature extraction algorithm for palm-print recognition[J].Digital Signal Processing,2013,23(1):244-258.
[2] TOUBIN M,DUMONT C,VERRECCHIA E P,et al.Multi-scale analysis of shell growth increments using wavelet transform[J].Computers & Geosciences,1999,25(8):877-885.
[3] 臧玉萍,張德江,王維正.小波分層閾值降噪法及其在發(fā)動機振動信號分析中的應用[J].振動與沖擊,2009,29(8):57-60.ZANG Yu-ping,ZHANG De-jiang,WANG Wei-zheng.Per-level threshold de-noising method using wavelet and its application in engine vibration analysis[J].Journal of Vibration and Shock,2009,29(8):57-60.
[4] 肖紅兵,楊錦舟,鞠曉東,等.V系統(tǒng)在隨鉆聲波測井數據降噪中的應用[J].中國石油大學學報:自然科學版,2009,33(2):58-62.XIAO Hong-bing,YANG Jin-zhou,JU Xiao-dong,et al.Application of V-system in acoustic logging while drilling data de-noising[J].Journal of China University of Petroleum(Edition of Natural Science),2009,33(2):58-62.
[5] ALONSO F J,DEL CASTILLO J M,PINTADO P.Application of singular spectrum analysis to the smoothing of raw kinematic signals[J].Journal of Biomechanics,2005,38(5):1085-1092.
[6] ZHAO X,YE B.Selection of effective singular values using difference spectrum and its application to fault diagnosis of headstock[J].Mech Syst Signal Process,2011,25(5):1617-1631.
[7] JIANG Y H,TANG B P,QIN Y,et al.Feature extraction method of wind turbine based on adaptive Morlet wavelet and SVD[J].Renewable Energy,2011,36(8):2146-2153.
[8] 段禮祥,張來斌,王朝暉,等.柴油機振動信號的小波包奇異值降噪[J].中國石油大學學報:自然科學版,2006,30(1):93-97.DUAN Li-xiang,ZHANG Lai-bin,WANG Zhao-hui,et al.De-noising of diesel vibration signal using wavelet packet and singular value decomposition[J].Journal of China University of Petroleum(Edition of Natural Science),2006,30(1):93-97.