徐彥凱，雙 凱

(中國(guó)石油大學(xué)地球物理與信息工程學(xué)院，北京 102249)

在油氣井下電磁脈沖數(shù)據(jù)傳輸中，有效信號(hào)具有暫態(tài)瞬變的特征。由于傳輸距離遠(yuǎn)，接收到的微弱有用信號(hào)往往被環(huán)境噪聲所淹沒。因此，有效檢測(cè)瞬變?nèi)跣盘?hào)是電磁脈沖數(shù)據(jù)傳輸?shù)年P(guān)鍵技術(shù)。小波分析是近些年發(fā)展起來的一種新的時(shí)頻分析方法，小波變換能使信號(hào)能量集中在一些大的小波系數(shù)中;而噪聲能量卻分布于整個(gè)小波域內(nèi)，對(duì)應(yīng)的小波系數(shù)小。因此，小波閾值法可以較好地用于抑制噪聲并檢測(cè)信號(hào)?？紤]到加噪信號(hào)對(duì)應(yīng)矩陣的奇異值分解一定程度上可以反映矩陣的主要特征，實(shí)現(xiàn)噪聲去除和信號(hào)提取，筆者將小波及奇異值分解相結(jié)合，實(shí)現(xiàn)微弱瞬變信號(hào)的檢測(cè)。

1 小波變換及降噪應(yīng)用

1.1 小波變換理論

設(shè)ψ(t)∈L2(R)，若ψ(t)為母小波，其傅里葉變換Ψ(ω)滿足條件:

將ψ(t)進(jìn)行伸縮和平移，得到小波基函數(shù)

其中a為尺度(伸縮)因子，b為平移因子。信號(hào)f(t)的連續(xù)小波變換定義［1］為

對(duì) a、b 進(jìn)行二進(jìn)離散化，即 a=2j，b=2jk，j，k∈Z，得到的二進(jìn)小波變換為

其逆變換為

Mallat算法實(shí)現(xiàn)了二進(jìn)小波的快速分解與重構(gòu)［2］，分解公式為

重構(gòu)公式為

式中，cj，k為尺度系數(shù);dj，k為小波系數(shù);h、g、～h、～g 分別為低通和高通分解與重建濾波器;j為分解層數(shù)。

1.2 小波閾值法降噪檢測(cè)技術(shù)

由于小波變換能使信號(hào)的能量集中在一些大的小波系數(shù)中，而噪聲能量卻平均分布于整個(gè)小波域內(nèi)，對(duì)應(yīng)的小波系數(shù)小。應(yīng)用信號(hào)與噪聲小波系數(shù)的特點(diǎn)對(duì)信號(hào)降噪，降噪步驟［3］如下:

(1)確定小波基函數(shù)和小波分解最大層數(shù)，應(yīng)用式(5)對(duì)信號(hào)進(jìn)行J層小波分解得到小波系數(shù)和尺度系數(shù)。

(2)確定每層小波系數(shù)的閾值，根據(jù)每層閾值按閾值函數(shù)形式對(duì)小波系數(shù)進(jìn)行處理，把保留下來的作為降噪后信號(hào)的小波系數(shù)。

(3)利用式(6)對(duì)信號(hào)進(jìn)行重構(gòu)，并對(duì)降噪后信號(hào)檢測(cè)。

上述步驟中，分解尺度、閾值函數(shù)和小波基函數(shù)的合理選擇對(duì)信號(hào)檢測(cè)至關(guān)重要。

2j尺度下小波偽頻率f2j與小波中心頻率fc的關(guān)系［3］:

式中，Δt為采樣周期。

若信號(hào)的中心頻率接近小波的偽頻率，則會(huì)使此頻率段的小波系數(shù)達(dá)到最大，從而有效地檢測(cè)信號(hào)。研究的瞬變信號(hào)的頻率集中在(0.3～1.1)×105Hz，db2小波的中心頻率為0.6667 Hz，采樣周期為1×10－6s，由式(7)可知，j=3和4的db2小波偽頻率在信號(hào)的頻率范圍，所以確定小波分解最大層數(shù)為4。

常用的方法有硬閾值函數(shù)和軟閾值函數(shù)，采用軟閾值處理會(huì)對(duì)所有的小波系數(shù)都進(jìn)行抑制，可使降噪后信號(hào)比較平滑，但對(duì)于提取具有突變特性的瞬變信號(hào)而言，對(duì)降噪后信號(hào)的平滑性沒有要求，相反采用硬閾值方法能夠更好地保留原信號(hào)的一些瞬變特征以便信號(hào)檢測(cè)，所以采用硬閾值處理方法更為合適。

選擇小波基函數(shù)時(shí)考慮正交小波變換后的數(shù)據(jù)不會(huì)出現(xiàn)冗余，并保證信號(hào)精確重構(gòu)。因此選擇幾種有代表性的常用正交小波基:haar小波、coif1小波、sym2小波和db2小波，在此基礎(chǔ)上進(jìn)一步通過仿真試驗(yàn)對(duì)比選優(yōu)。

1.3 仿真與分析

四種小波基下的重建信號(hào)見圖1。

圖1 四種小波基下的重建信號(hào)Fig.1 Reconstructed waveforms by four wavelet functiones

電磁脈沖傳輸?shù)臄?shù)字信號(hào)用圖1(a)中的瞬變信號(hào)表示，該信號(hào)是一個(gè)持續(xù)時(shí)間短的暫態(tài)過程。通過電路分析，得到該信號(hào)的數(shù)學(xué)模型為

式中，A為信號(hào)幅度;ω為振蕩頻率;α為衰減因子，決定信號(hào)暫態(tài)過程的持續(xù)時(shí)間，其持續(xù)時(shí)間由電路參數(shù)決定;t≤t0為充電時(shí)間;t＞t0為放電時(shí)間。α=15，ω =4.7 × 105rad/s。

井下電磁脈沖信號(hào)傳輸中，噪聲源主要有兩大類:井下鉆井環(huán)境產(chǎn)生的噪聲和電子器件產(chǎn)生的電子噪聲。鉆井噪聲頻率范圍為1～4 kHz。電子噪聲在頻域和時(shí)域上分布一致，是一種高斯白噪聲［4］。由于鉆井噪聲可通過濾波的方法去噪，在此只考慮無法直接濾除的白噪聲。

針對(duì)小波基的選取進(jìn)行試驗(yàn)，試驗(yàn)中采用如圖1(b)所示的－2.5 dB加噪信號(hào)，統(tǒng)一的4尺度分解，得到圖1所示的4種小波基下的重構(gòu)信號(hào)。比較圖1中的(c)、(d)、(e)、(f)，可以直觀地看出，針對(duì)本文研究的瞬變信號(hào)，db2小波的降噪效果最佳，sym2次之，coif1小波最差，因此選擇db2小波基。

由圖1(b)信號(hào)作4尺度db2小波分解的各尺度系數(shù)如圖2所示。從圖2看出，無法從j=1，2的小波系數(shù)中區(qū)分信號(hào)和噪聲，因此重構(gòu)時(shí)將j=1，2的小波系數(shù)全置零，保留 j=3，4中大于閾值 λ=σ的小波系數(shù)。

圖2 db2小波分解各尺度系數(shù)Fig.2 Scale coefficients with db2 wavelet

圖1(c)表明，在高信噪比的情況下，小波閾值法可以實(shí)現(xiàn)信號(hào)的檢測(cè)，但從圖3所示的 －8.5 dB加噪信號(hào)的db2小波重建信號(hào)可以看出:在(0.1～0.2)ms和(0.7～0.8)ms處沒有信號(hào)，但在相應(yīng)位置重構(gòu)信號(hào)的幅度較大。

因此，為了更好地檢測(cè)強(qiáng)噪聲中的弱信號(hào)，需要將j=3，4信號(hào)和噪聲的小波系數(shù)進(jìn)一步分離。

圖3 -8.5 dB加噪信號(hào)的db2小波重建信號(hào)Fig.3 db2 WT reconstructed waveform with-8.5 dB SNR

2 奇異值分解降噪

2.1 奇異值分解理論

設(shè)矩陣Hm×n(m≤n)秩為r，則存在兩個(gè)正交矩陣U、V和對(duì)角矩陣D，使得下式成立［8］:

式中，λi為矩陣H的奇異值;ui和vi是λi的特征向量為特征項(xiàng)。奇異值分解是一種代數(shù)特征提取方法，能找到矩陣各向量之間的內(nèi)在本質(zhì)屬性。因此加噪信號(hào)序列的奇異值分解可以提取序列主要特征，實(shí)現(xiàn)信號(hào)降噪和檢測(cè)。其步驟如下:

(1)Hankle矩陣生成。設(shè)加噪信號(hào)序列為X=［x0，x1，…，xN－1］，構(gòu)造 Hankle 矩陣 H［5-7］

矩陣H的行數(shù)L、列數(shù)J及序列X的長(zhǎng)度N之間應(yīng)滿足

(2)矩陣奇異值分解。對(duì)矩陣H進(jìn)行奇異值分解，λ1，λ2，…，λr為矩陣 H 的按降序排列的非零奇異值，即 λ1≥ λ2≥ … ≥ λr。

由式(12)重建信號(hào)并檢測(cè)有用信號(hào)存在與否。

2.2 參數(shù)選擇

奇異值分解中，窗長(zhǎng)L及保留特征值個(gè)數(shù)TR的選擇影響降噪和檢測(cè)效果。

本文首先選取TR=5，窗長(zhǎng)L分別選取15、22、30、40、60，碼元中有用信號(hào)的采樣點(diǎn)數(shù)是 30，加噪信號(hào)的信噪比分別為－2.5、－5.5、－8.5、－9.5、－10.5、－11.5 dB，仿真得到如圖4所示的SVD降噪峰值信噪比與窗長(zhǎng)L的關(guān)系。

圖4 窗長(zhǎng)對(duì)SVD降噪峰值信噪比的影響Fig.4 Effect of window length on SVD-denoised PSNR

從圖4看出:當(dāng)原信噪比大于－8.5 dB，窗長(zhǎng)L選取30效果最好;當(dāng)原信噪比小于－8.5 dB，窗長(zhǎng)L選取22效果最好;窗長(zhǎng)L選取過大(L=60)或過小(L=15)降噪效果明顯差。因此得出:①若原信噪比較大，選擇窗長(zhǎng)L等于碼元中有用信號(hào)長(zhǎng)度時(shí)效果最佳;②若原信噪比較小，選擇窗長(zhǎng)L比碼元中有用信號(hào)的長(zhǎng)度稍小些降噪和檢測(cè)效果最好，這是由于信號(hào)末端采樣值接近零，當(dāng)噪聲較強(qiáng)時(shí)，奇異值分解中這部分值的特征完全被噪聲淹沒造成的。

然后選取窗長(zhǎng)L=22，保留特征值個(gè)數(shù)TR分別選取15、10、5、3、2，加噪信號(hào)的信噪比分別為 －2.5、－5.5、－8.5、－9.5、－10.5、－11.5 dB，得到如圖5所示的SVD降噪峰值信噪比與保留特征值個(gè)數(shù)TR的關(guān)系。

圖5 保留特征值個(gè)數(shù)TR對(duì)SVD降噪峰值信噪比的影響Fig.5 Effect of main singular values1 number on SVD-denoised PSNR

從圖5看出:①降噪后峰值信噪比隨著原信噪比的減小而減小;②保留特征值個(gè)數(shù)TR=2峰值信噪比最高，檢測(cè)效果也最佳;③原信噪比越小，TR的選擇對(duì)降噪后峰值信噪比的影響越大，這是由于奇異值分解中大部分特征值已經(jīng)體現(xiàn)不出信號(hào)的特性。

2.3 仿真與分析

在研究奇異值分解參數(shù)選擇的基礎(chǔ)上，采用奇異值分解方法對(duì)－8.5 dB加噪信號(hào)檢測(cè)，其中選擇參數(shù)L=22，TR=2，仿真結(jié)果如圖6所示。

圖6 -8.5 dB加噪信號(hào)的SVD重建信號(hào)Fig.6 Reconstructed waveform by SVD method with-8.5 dB SNR

比較圖3(c)和圖6(c)可知:從檢測(cè)的角度看，奇異值分解比小波降噪后判決出錯(cuò)的概率小;但從輸出信噪比看，奇異值分解比小波降噪差。

3 小波奇異值分解降噪檢測(cè)

通過上述分析可知，若將小波和奇異值相結(jié)合［7-8］會(huì)使信號(hào)和噪聲的小波系數(shù)更好地分離，檢測(cè)效果更佳。具體方法是:首先將信號(hào)作小波分解，再對(duì)小波分解系數(shù)作奇異值分解，最后通過閾值法保留信號(hào)小波系數(shù)并重建降噪信號(hào)，利用重建信號(hào)進(jìn)行信號(hào)檢測(cè)。

圖7(a)和圖7(b)分別為信號(hào)和－8.5 dB加噪信號(hào)，采用小波閾值、奇異值分解及小波奇異值三種方法對(duì)圖7(b)加噪信號(hào)降噪，降噪結(jié)果分別如圖7(c)～(e)?？梢钥闯?僅采用奇異值分解的圖7(d)降噪效果最差;采用小波降噪的圖7(c)雖然效果較好，但在沒有信號(hào)的(0.1～0.4)ms、和(0.5～0.6)ms處的重構(gòu)信號(hào)的幅度較大;采用小波奇異值方法的圖7(e)則可以準(zhǔn)確檢測(cè)信號(hào)的存在。因此與小波閾值法和奇異值法相比，小波奇異值法降噪和檢測(cè)效果最佳。

圖7 -8.5 dB加噪信號(hào)的三種方法重建信號(hào)Fig.7 Reconstructed waveforms by three methods with-8.5 dB SNR

4 結(jié)束語

針對(duì)油氣井下電磁脈沖數(shù)據(jù)傳輸中瞬變?nèi)跣盘?hào)的檢測(cè)，研究了小波閾值法和奇異值分解及參數(shù)選擇。在此基礎(chǔ)上提出了小波與奇異值分解相結(jié)合降噪檢測(cè)信號(hào)的方法。仿真結(jié)果表明，該方法能更好地區(qū)分信號(hào)和噪聲，獲得了更好的降噪和檢測(cè)結(jié)果。

［1］ IMTIAZ H，F(xiàn)ATTAH 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ā)動(dòng)機(jī)振動(dòng)信號(hào)分析中的應(yīng)用［J］.振動(dòng)與沖擊，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)在隨鉆聲波測(cè)井?dāng)?shù)據(jù)降噪中的應(yīng)用［J］.中國(guó)石油大學(xué)學(xué)報(bào):自然科學(xué)版，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ī)振動(dòng)信號(hào)的小波包奇異值降噪［J］.中國(guó)石油大學(xué)學(xué)報(bào):自然科學(xué)版，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.